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Abstract 


The  protein  folding  problem  attempts  to  predict  the  native,  or  folded,  state  of  a 
protein  in  three-dimensional  space,  given  its  primary  sequence  of  amino  acids.  One  com¬ 
mon  approach  for  a  solution  is  to  treat  each  complex  amino  acid  as  a  single  sphere,  or 
“united  atom,”  and  to  model  each  peptide  linkage  between  residues  by  a  virtual  bond 
between  spheres.  Computational  efforts  being  examined  rely  on  two  major  assumptions: 
1)  for  any  specific  molecular  conformation,  a  corresponding  potential  energy  function  can 
be  computed,  and  2)  the  three-dimensional,  folded  state  corresponds  to  the  global  mini¬ 
mum  of  this  energy  function.  The  optimization  method  being  used  to  minimize  the 
potential  energy  involves  collecting  a  large  number  of  conformers,  each  attained  by  find¬ 
ing  a  local  minimum  of  the  potential  energy  function  from  a  random  starting  point.  The 
information  from  these  conformers  is  then  used  to  form  a  convex  quadratic  global  under¬ 
estimating  function  for  the  potential  energy  of  the  known  conformers.  The  minimum  of 
this  underestimator  is  used  to  predict  the  global  minimum  for  the  function,  allowing  a 
localized  conformer  search  to  be  performed  based  on  the  predicted  minimum.  The  new  set 
of  conformers  generated  by  the  localized  search  can  serve  as  the  basis  for  another  qua¬ 
dratic  underestimation.  After  several  repetitions,  the  global  minimum  can  be  found  with 
reasonable  assurance.  The  conformer  which  lies  at  the  global  minimum  represents  the 
three-dimensional  folded  state  of  the  molecule. 
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I.  Introduction 

A.  The  Protein  Folding  Problem 

One  of  the  most  important  and  most  difficult  problems  in  biophysics  and  biochem¬ 
istry  is  the  protein  folding  problem.  In  the  late  1950’s  Christian  B.  Anfinsen  and  his  col¬ 
leagues  at  the  National  Institute  of  Health  made  the  remarkable  discovery  that  the  folded 
state  of  a  protein  was  completely  dependent  on  the  one-dimensional  linear  sequence  (i.  e. 
“primary”  sequence)  of  amino  acids  from  which  the  protein  is  constructed.  External  fac¬ 
tors,  such  as  enzymes  present  at  the  time  of  folding,  have  no  effect  on  the  final  or  native 
state  of  the  protein.  This  led  to  the  formulation  of  the  protein  folding  problem:  Given  a 
known  primary  sequence  of  amino  acids,  what  would  be  its  native,  or  folded,  state  in  3- 
dimensional  space?  Despite  over  thirty  years  of  research,  a  solution  to  the  protein  folding 
problem  remains  elusive  (Richards  1991). 

A  solution  to  the  protein  folding  problem  would  provide  immeasurable  benefits  to 
the  biotechnology  field.  Many  new  products  in  the  biotechnology  industry  are  novel  pro¬ 
teins.  It  is  already  possible  to  design  “custom”  amino  acid  sequences  quickly  with  DNA 
sequence  analysis;  however,  the  inability  to  predict  the  folded  state  of  the  protein,  which 
determines  the  protein’s  functionality,  forces  biotechnology  industries  to  rely  on  long, 
expensive  methods  of  development.  When  Eli  Lilly  made  his  first  attempts  at  producing 
insulin  for  diabetes,  a  biologically  inactive  material  resulted,  which  resembled  scrambled 
eggs.  While  he  was  able  to  splice  the  human  gene  for  insulin  into  DNA  molecules  of  bac¬ 
teria  and  get  them  to  synthesize  the  needed  protein  chain,  incorrect  folding  caused  the  pro¬ 
cess  to  fail.  Only  after  considerable  expenses  in  time  and  money  was  the  correct  folded 
sequence  found.  Currently,  the  three-dimensional  folded  state  of  a  protein  can  only  be 
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found  using  X-ray  crystallography,  an  enormously  expensive  and  time  consuming  process 
which  can  take  as  much  as  several  months  or  even  years. 

While  methods  have  been  developed  which  attempt  to  predict  a  folded  protein’s 
secondary  structure  (distinct  sections  having  characteristic  shapes  including  helices,  beta 
strands,  and  turns),  these  methods  rarely  yield  information  which  provides  a  description  of 
the  tertiary  structure  (the  folded  state).  As  a  result,  the  protein  folding  problem  has  gener¬ 
ally  been  thought  to  be,  if  not  completely  insoluble,  not  likely  to  be  solved  for  a  very  long 
time. 

B.  Possible  Solutions  to  the  Protein  Folding  Problem 

Recent  events  have  renewed  hope  for  a  solution  to  the  problem.  Several  successful 

predictions  have  been  made  and  announced  before  the  experimental  structures  were 
known.  While  most  of  these  predictions  have  been  made  with  a  blend  of  a  human  expert’s 
abilities  and  computer  assistance,  fully  automated  methods  have  shown  promise  for  pro¬ 
ducing  previously  unattainable  accuracy.  Advances  in  theory,  experiment,  and  computing 
power,  coupled  with  growing  interest  from  industry,  have  kindled  new  optimism  for  a 
solution  to  the  protein  folding  problem  (Richards  1991,  Benner  and  Gerloff  1993). 

Computational  attempts  to  solve  the  folding  problem  can  be  divided  into  two  basic 
strategies:  human  based  and  machine  based  predictions.  The  first  relies  on  a  close  interac¬ 
tion  between  a  human  expert  and  a  computer  model  which  simulates  the  behavior  of  the 
protein.  One  such  system  was  developed  by  Mark  C.  Surles  at  the  University  of  North 
Carolina  at  Chapel  Hill.  The  graphical  modeling  system  he  developed,  “Sculpt”,  pro¬ 
vides  an  environment  in  which  an  expert  can  graphically  manipulate  the  protein  into  a 
folded  state.  Through  the  use  of  the  Sculpt  program,  the  expert  can  insure  that  the  protein 
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structure  complies  with  established  biochemical  rules  (Surles  1992).  While  this  approach 
is  a  huge  improvement  over  a  purely  human  based  analysis  of  protein  structure,  it  still 
remains  dependent  on  possibly  enormous  amounts  of  time  spent  by  an  expert. 

Machine  based  prediction  strategies  attempt  to  lessen  the  reliance  on  experts  by 
developing  a  completely  computational  method.  Such  approaches  are  generally  based  on 
two  assumptions.  First,  there  exists  a  potential  energy  function  for  the  protein;  and  second, 
the  folded  state  corresponds  to  the  structure  with  the  lowest  potential  energy  (minimum  of 
the  potential  energy  function)  and  is  thus  in  a  state  of  thermodynamic  equilibrium.  This 
view  is  supported  by  in-vitro  observations  that  proteins  can  successfully  refold  from  a 
variety  of  denatured  states.  The  alternate  assumption  is  that  the  native  state  of  the  protein 
may  be  at  a  local  energy  minimum  which  can  be  separated  from  the  true  global  energy 
minimum  by  a  large  kinetic  barrier.  Recent  studies  of  the  a-Lytic  protease  revealed  two 
conformational  forms,  active  and  inactive,  separated  by  just  such  a  kinetic  barrier  (Baker 
et  al.  1992).  While  the  active,  or  stable,  form  of  the  molecule  occurred  at  the  global  mini¬ 
mum,  a  catalyst  is  required  to  change  the  “molten  globule”  found  at  the  inactive  state  into 
the  final  state  for  the  protein.  This  evidence  seems  to  support  the  view  that  the  folded  state 
may  not  lie  at  the  global  minimum  in  every  case.  Ruben  A.  Abagyan  raises  two  questions 
related  to  this  phenomenon  (Abagyan  1993): 

( 1 )  Could  a  metastable  state  have  all  the  features  of  the  normal  protein? 

(2)  If  it  could,  is  this  behavior  typical  for  all  proteins? 

Currently  there  is  no  clear  evidence  for  a  positive  answer  to  the  first  question.  The  a-Lytic 
experiments,  while  raising  important  questions,  do  not  provide  evidence  concerning  the 
properties  of  the  normal  protein  (the  final  state  for  a-Lytic  occurred  at  the  global  mini- 
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mum).  The  second  question  can  be  answered  definitively  in  the  negative.  Small  single¬ 
domain  proteins  do  fall  into  their  lowest  energy  state.  Evolutionary  theory  also  supports  a 
folded  state  at  a  global  energy  minimum.  Protein  sequences  have  evolved  under  pressure 
to  perform  certain  functions,  which  for  most  known  occurrences  require  stable,  unique 
and  compact  structures.  Unless  specifically  required  for  a  certain  function,  there  was  no 
biochemical  need  for  proteins  to  hide  their  global  minimum  behind  a  large  kinetic  energy 
barrier.  While  kinetic  blocks  may  occur,  they  should  be  limited  to  special  proteins  devel¬ 
oped  for  certain  functions,  or  occur  randomly  with  the  probability  of  kinetic  blocks  grow¬ 
ing  with  problem  size  (Abagyan  1993). 

C.  Models  for  Computational  Methods 

Unfortunately,  finding  the  “true”  energy  function  of  a  protein,  if  one  even  exists,  is 

so  unimaginably  complicated  that  there  are  no  means  nor  hope  of  determining  it  exactly. 
With  proteins  ranging  in  size  up  to  1,053  amino  acids  (a  collagen  found  in  tendons),  and 
covering  an  almost  infinite  number  of  possible  conformations,  exhaustive  conformational 
searches  will  never  be  possible. 

One  possible  way  of  finding  the  global  energy  minimum  is  to  use  a  simplified 
model  of  the  protein  system.  By  using  a  simplified  model,  the  complexity  of  the  problem 
formulation  could  be  reduced  to  an  acceptable  level  for  optimization  techniques.  Care 
must  be  taken,  however,  to  insure  that  the  error  included  in  such  an  approximation  does 
not  drive  the  computational  solution  too  far  away  from  the  true  native  state.  Simplified 
computational  protein  models  can  be  currently  divided  into  two  major  classes:  lattice 
restricted  models  and  models  in  continuous  space. 

The  first  of  these  methods,  lattice  models,  relies  on  constraining  the  position  of  the 


10 


'mino  acids  to  the  vertices  of  a  three-dimensional  lattice.  By  limiting  the  possible  posi¬ 
tions  for  amino  acid  placement,  lattice  methods  have  the  potential  for  large  increases  in 
algorithm  speed.  This  increase,  however,  is  reliant  on  the  quality  of  the  grid  developed.  A 
lattice  with  a  very  large  number  of  vertices  could  provide  an  accurate  model  of  the  protein. 
Unfortunately,  this  increase  in  accuracy  carries  a  corresponding  decrease  in  speed.  Con¬ 
versely,  a  very  coarse  lattice  with  relatively  few  vertices  enables  fast  execution  time,  but 
may  not  provide  an  accurate  solution. 

The  second  class  of  techniques  attempts  to  avoid  the  reliance  on  lattice  restrictions 
by  performing  the  optimization  in  continuous  space.  These  methods  represent  a  protein  as 
a  string  of  beads  (amino  acids)  whose  position  is  defined  by  their  location  relative  to  the 
previous  three  beads.  One  significant  advantage  of  this  formulation  of  the  folding  problem 
is  that  it  allows  the  model  to  take  advantage  of  known  scientific  knowledge  about  the 
chemical  construction  of  proteins.  The  use  of  knowledge  such  as  the  Ramachandran  plot 
(Lehninger  1970)  which  specifies  the  possible  values  for  certain  angles  in  a  protein,  fur¬ 
ther  simplifies  the  problem.  This  class  of  problems  is  also  appealing  since  the  problems 
may  be  subjected  to  a  variety  of  global  optimization  methods  including  minimization  from 
random  starting  conformations,  molecular  dynamics,  Monte  Carlo  techniques,  and  simu¬ 
lated  annealing  (Abagyan  1993). 

This  paper  will  examine  several  possibilities  for  the  global  optimization  of  a  pro¬ 
tein’s  potential  energy  function.  The  approaches  rely  on  two  major  assumptions: 

1 .  For  any  specific  molecular  conformation,  a  corresponding  potential  energy 
can  be  computed  from  a  potential  energy  function  model  of  the  protein. 

2.  The  native  state  corresponds  to  the  global  (or  near  global)  minimum  of  this 
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energy  function. 

Each  amino  acid  is  represented  as  a  single  point  whose  position  is  determined  by  a  bond 
length,  bond  angle,  and  torsion  angle  triple.  Several  different  methods  for  this  optimiza¬ 
tion,  including  both  continuous  and  lattice  based  approaches,  will  be  examined. 
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II.  The  Potential  Energy  Function 

A.  Molecular  Representations 

Molecular  structure  information  is  generally  given  in  terms  of  bond  lengths,  bond 
angles,  and  torsion  angles.  A  protein  consists  of  n  amino  acids  in  the  “primary”  sequence, 
aj,a2,...an,  where  a,  represents  the  Ith  amino  acid  in  the  sequence.  For  every  consecutive 
pair  of  amino  acids,  a,-.]  and  a,,  /,  represents  the  distance  by  which  they  are  separated.  For 
every  three  consecutive  amino  acids,  a,-_2,  a,_i,  a,,  the  bond  angle,  0(,  represents  the  posi¬ 
tion  of  the  third  amino  acid  with  respect  to  the  line  containing  the  previous  two  amino 
acids.  Similarly,  for  every  four  consecutive  amino  acids,  a,_3,  a(.2,  a,.j,  and  a,,  the  torsion 
angle,  cp(,  represents  the  relative  position  of  the  fourth  amino  acid,  a,,  with  respect  to  the 
plane  containing  the  previous  three  amino  acids  (See  Figure  2.1). 


Figure  2.1  Conformation  for  a  Four  Amino  Acid 
Sequence 


B.  Potential  Energy  Function 

Empirical  representations  of  the  protein’s  potential  energy  include  energy  terms  to 
represent  chemical  bonds,  angles,  and  torsions,  as  well  as  non-bonded  interactions 
between  amino  acids  farther  apart  in  the  chemical  structure.  This  simplified  definition 
resembles  a  physical  model  in  which  beads  (amino  acids)  are  connected  by  springs  (first 
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term  of  Equation  2.1)  at  a  distance  of  /,.  Bond  angles,  0,,  which  are  determined  by  a 
sequence  of  three  amino  acids  (a^.a^.a,)  are  maintained  by  similar  “springs”  (the  second 
term).  Torsion  angles,  <p(,  are  modeled  by  a  trigonometric  based  penalty  (third  term  of 
Equation  2.1).  Such  a  formulation  is  often  called  a  molecular  mechanics  potential.  There¬ 
fore,  one  formulation  frequently  shown  to  express  the  overall  potential  is  shown  in  Equa¬ 
tion  2.1,  where  k[  and  are  the  bond  stretching  and  angle  bending  force  constants 
respectively.  In  the  third  term  of  the  equation,  V~  is  the  n-fold  torsional  constant  with  a 
phase  shift  of  8.  This  term  provides  for  preferred  torsion  angles  for  the  protein.  For  exam¬ 
ple,  if  n  =  3  and  the  phase  shift,  5,  is  zero,  preferred  torsion  angles  (configurations  result¬ 
ing  in  a  small  penalty)  can  be  found  at  60°,  180°,  and  300°.  (Since  the  cosine  of  180°  or 
any  odd  multiple  of  180°  is  -1,  the  penalty  incurred  by  the  preferred  torsions  will  be  zero.) 
The  constants,  1$  and  0O  represent  the  preferred  bond  length  and  bond  angle  for  each  set  of 
amino  acids. 


U  =  X^/^i-/o)2  +  S^e^0«~  ®o)2  +  Sy  l1  +  cos(n<p/  +  5)J  Equation  2. 1 
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The  final  portion  of  the  equation,  ^e.  |  —  |  -2j  —  |  |,  is  known  as  the 

«,  / 

Lennard-Jones  pairwise  potential.  This  term  defines  the  potential  energy  contributions  of 
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amino  acids  separated  by  more  than  two  residues  along  the  primary  chain.  Graphically,  the 
Lennard-Jones  term  can  be  represented  as  a  plot  (Figure  2.2)  of  the  potential  energy  of  the 
protein  (in  kcal)  versus  the  distance  (in  A)  separating  the  amino  acid  residues  in  the  pro- 


14 


tein.  In  Equation  2.1,  e,-  /  and  a,  /  are  the  Lennard- Jones  coefficients,  which  are  constants 


Figure  2.2  Lennard- Jones  Pairwise  Potential  (£=0.181,  a=4.0) 
defined  by  the  relationships  between  the  two  specific  beads  (amino  acids)  involved.  The 
r,  i  terms  in  the  Lennard-Jones  expression  represent  the  Euclidean  distance  between  the 
amino  acids  a,  and  at  (Ferguson  et  al.  1994)  (Case  1993). 


III.  Conversions 


A.  Coordinate  Representations 

Potential  energy  functions  for  molecular  conformations  usually  involve  computing 

the  distance  between  molecules.  For  example,  the  potential  function  discussed  in  Chapter 

((C.  ,V2 

II  (page  12)  contains  the  Lennard- Jones  pairwise  potential,  £e.  t  —  -2 

i,l  ’  \Vri,l 

which  uses  the  distance  between  all  amino  acids  separated  by  more  than  two  residues 
along  the  primary  chain  (such  as  amino  acids  aj  and  a4).  Unfortunately,  computing  dis¬ 
tances  using  bond  lengths,  bond  angles,  and  torsion  angles  is  extremely  difficult  and  very 
time  consuming.  Since  optimization  requires  this  computation  to  be  executed  often,  it  is 
much  more  profitable  to  convert  the  problem  into  cartesian  coordinates  before  performing 
the  optimization.  Once  this  conversion  has  been  completed,  the  length  between  the  resi¬ 
dues  can  be  found  with  a  minimal  amount  of  computation  time.  However,  the  conversion 
itself  becomes  very  complex  since  every  set  of  three  residues  creates  its  own  reference 
frame  for  bond  lengths,  bond  angles,  and  torsion  angles. 

B.  The  Conversion  Geometry 

To  perform  the  conversion  from  a  chemical  representation  (bond  lengths,  bond 
angles,  and  torsion  angles)  to  a  cartesian  representation,  the  first  three  points  in  the 
sequence  can  be  fixed,  without  loss  of  generality,  as  depicted  in  Figure  2. 1 .  The  first  amino 
acid,  a],  is  fixed  at  the  origin,  (0,0,0).  The  second  amino  acid,  a2,  is  positioned  at  (/2,0,0), 
which  is  a  distance  from  aj  equal  to  the  bond  length  (/2)  along  the  x  axis.  The  location  of 
the  third  amino  acid,  a3,  is  fixed  at  (/2+/3Cos(03),/3sin(03),O).  The  position  of  the  fourth 
and  any  other  amino  acids  constituting  the  primary  sequence  are  found  using  the  cartesian 
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»4 


Figure  3.1  Four  Amino  Acid  Sequence  (Cartesian  Coordinates) 

data  associated  with  the  previous  three  points,  and  the  bond  length,  bond  angle,  and  tor¬ 
sion  angle  associated  with  the  amino  acid  under  consideration. 

The  position  of  a  fourth  point,  a/,  in  a  four  amino  acid  sequence,  a;,  a ;,  ak,  a/,  is 
completely  determined  by  its  chemical  representation:  bond  length  (//),  bond  angle  (0/), 
and  torsion  angle  (<fy).  To  convert  this  information  into  x,  y,  and  z  coordinates,  the  follow¬ 
ing  quantities  must  be  computed:  (1)  a  plane  based  on  the  bond  angle  and  bond  length,  (2) 
a  plane  based  on  the  torsion  angle  and  bond  length,  and  (3)  a  sphere  with  radius  equal  to 
the  bond  length.  The  intersection  of  the  two  planes  yields  a  line,  which  can  then  be  inter¬ 
sected  with  the  sphere  to  provide  two  possible  positions  for  the  bead  a/.  The  ambiguity 
between  these  two  points  can  be  resolved  by  closer  inspection  of  the  definition  of  the  tor¬ 
sion  angle.  The  remaining  sections  of  this  chapter  serve  to  illustrate  this  technique. 

C.  Plane  based  on  Bond  Angle 

The  bond  angle,  0/,  and  bond  length,  //,  define  a  plane  in  which  the  fourth  amino 
acid,  a/,  must  lie  (See  Figure  3.2).  The  normal  vector  to  this  plane  is  defined  as  the  vector 
between  the  two  previous  amino  acids  (in  Figure  3.2,  a j  and  ak).  If  rk  represents  the  posi¬ 
tion  of  atom  k,  and  r,  represents  the  position  of  atom  j,  then: 
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Figure  3.2  Plane  Found  From  Bond  Angle  and  Bond  Length 


r kj  =  r k-Tj  Equation  3.1 

and  this  vector,  in  fact,  is  the  normal  vector,  nQ,  to  the  plane,  Pe.  In  addition  to  a  normal 
vector,  a  point  on  the  plane  must  be  specified  to  complete  the  construction  of  the  plane. 
The  bond  angle,  when  combined  with  the  bond  length,  provides  enough  information  to 
compute  this  point.  That  is,  the  distance  from  the  previous  point,  a*,  to  the  plane  is  given 
by: 

d  =  -/, cosQr  Equation  3.2 

Since  this  distance  is  taken  in  the  direction  of  the  normal  vector,  a  point  on  the  plane  is 
given  by: 

ne 

rQ  =  d n — ii  +  r, .  Equation  3.3 

INI 

Therefore,  the  equation  of  the  plane  is: 


Aex  +  fley  +  Cez  -  D0  =  0 


Equation  3.4 
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where  Ae,  Be,  and  Ce  are  the  i,  j,  and  k  components  of  the  normal  vector,  iiq.  The  constant 
term,  D0,  is  given  by: 

D0  =  Abxq  +  fl0ye  +  CqZq  Equation  3.5 

where  (xe,  ye,  ze)  is  the  point  re. 

D.  Plane  based  on  Torsion  Angle 

The  position  of  the  amino  acid,  a/,  can  be  further  restricted  by  a  plane,  P^,  con¬ 
structed  from  the  torsion  angle  and  bond  length  (see  Figure  3.4).  The  normal  vector  to  the 
plane,  P<p  is  found  by  taking  the  cross  product  of  rkj  and  (for  r kj  see  Equation  3.1)  so 


Figure  3.3  Plane  Derived  From  Torsion  Angle 
and  Bond  Length 


that: 


ru  =  r.-r, 


%  =  W 


Equation  3.6 


Equation  3.7 
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This  plane  is  will  be  orthogonal  to  the  first  plane,  since  its  normal  vector  implicitly 
involves  the  cross  product  of  the  bond  angle  plane’s  normal  vector  (ne  =  ).  A  point  on 

the  plane,  P^,  can  be  found  using  the  distance  from  a*  to  the  plane,  Pq,,  and  the  known  nor¬ 
mal  vector,  n^.  The  distance  to  the  plane  is  obtained  by: 

d  =  ^sin  ( 180°  -  0,)  sin  ((p;) .  Equation  3.8 


The  point,  iy  on  this  plane  is  then: 


Equation  3.9 


As  in  section  C,  the  equation  for  the  plane  is  given  by: 

+  Bfy  +  C^z  -  £><p  =  0  Equation  3.10 

where  A^,  Bq,,  and  Cq,  are  the  i,  j,  and  k  components  of  the  normal  vector,  liq,.  The  constant 
term,  Dq,,  is  given  by: 

D(P  =  A9*<p  +  V<p  +  C«PZ(P  Equation  3.11 

where  (Xq,,  y9,  z(?)  is  the  point  r«p. 

E.  Bond  Length  Sphere 

Finally,  a  sphere  of  radius  equal  to  the  bond  length  about  the  last  point,  a*,  is 
described.  The  equation  is  simply: 

(x-xk)2+  ( y-yk)2+  (z  -  zk) 2  =  If  Equation  3.12 


where  (x*,y*,z*)  is  the  cartesian  position  of  the  point  a*. 
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F.  Line  at  Intersection  of  Bond  Angle  and  Torsion  Angle  Planes 

Since  the  bond  angle  plane  and  torsion  angle  plane  are  orthogonal  and  therefore 

guaranteed  to  intersect,  the  line  formed  by  this  intersection  can  be  found.  The  direction  of 
this  line,  its  line  vector,  is: 

l  =  itgXn^.  Equation  3.13 


The  cartesian  components  of  this  line  vector  are  denoted  by,  (A/,  B/,  C /).  To  complete  the 
line  definition,  a  point  on  the  line  must  be  known.  This  can  be  easily  accomplished  by 
finding  any  solution  to  the  system  of  two  equations  defining  the  two  planes,  Pe  and  P^: 


Q 

A  B  r 

V 

D 

(p  (p  <p 

L  m 

Equation  3.14 


By  setting  any  single  variable  to  zero  (since  the  system  is  underdetermined)  a  point, 
rpt  =  (*p/»  ypr  zpt)  which  lies  on  the  line  can  be  obtained.  Therefore,  using  a  parametric 
formulation  (with  t  representing  the  parameter),  the  line  is: 

x  =  xpl  +  Att  Equation  3.15 

y  =  y  +  B(t  Equation  3.16 

z  =  zpt+Clt.  Equation  3.17 


G.  Intersection  of  Line  and  Sphere 

By  combining  the  information  from  the  line  found  in  the  previous  section  and  the 

sphere  based  on  bond  length,  the  position  of  the  new  residue  can  be  limited  to  two  possi¬ 
ble  points  (Figure  3.4).  Substituting  the  line  equations  into  the  sphere  equation  developed 
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in  Equation  3.12,  one  can  derive  the  following  relationships: 

(Xpt+Alt-xk)2+  (ypi  +  Blt-yk)\  (Zpt  +  C,t-zk)2  =  l) . 

Since  t  is  the  only  unknown,  by  the  quadratic  rule, 

-b±  *]b2  -4ac 
2  a 


where. 


a  =  +  B]  +  c]  J 


b  =  2[A,{xpt-xk)  +Bl(ypt-yk)  +C,(zpf-zt)] 
c=  (xP,-xk)2+  (ypt-yf+^p,-^2- 


Equation  3.18 

Equation  3.19 

Equation  3.20 
Equation  3.21 

Equation  3.22 


The  two  possible  points  can  be  easily  found  by  solving  Equation  3.15  through  Equation 
3. 17  using  the  two  values  of  t  derived  above. 
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H.  Resolving  to  One  Point 

The  ambiguity  for  the  two  points  can  be  resolved  by  simply  testing  the  resultant 
torsion  angles  against  the  desired  torsion  angle,  tpj.  One  of  the  two  points  is  selected  and 
its  torsion  angle  is  computed.  If  this  torsion  angle  is  within  a  small  error  of  the  desired  tor¬ 
sion  angle,  (p/,  the  point  is  accepted.  Otherwise,  the  point  is  discarded  and  the  other  point 
defined  by  the  intersection  is  selected.  Hence,  to  complete  the  computation  of  the  torsion 
angle  three  vectors,  (Equation  3.6),  rkj  (Equation  3.1),  and  r^,  are  defined: 

r[k  =  rl-rk.  Equation  3.23 

Two  new  vectors  are  constructed:  t,  which  is  normal  to  the  plane  defined  by  a,,  a /,  and  a*, 
and  u,  which  is  normal  to  the  plane  defined  by  a,,  a*,  a/: 

t  =  r  x  r .  Equation  3.24 

IJ  KJ 

u  =  rlk  x  rkj .  Equation  3.25 

The  torsion  angle,  (p,  is  defined  as  the  angle  between  the  two  planes  passing  through  a,,  a j, 
and  a*,  and  a /,  a*,  a/,  respectively.  Therefore: 

*  ■  aC°S(l«)-  Equation  3.26 

When  t  and  u  are  parallel,  <p  =  0.  If  t  and  u  are  antiparallel,  then  <p  =  180°  (Swope  and  Fer¬ 
guson  1992).  This  value  of  tp  is  compared,  as  described  above,  to  <p/  to  eliminate  any  ambi¬ 
guity  and  therefore  conclude  the  computation  of  the  cartesian  representation  of  the  point  a/ 
(for  another  method  of  computing  the  conversions,  see  Maranas  and  Floudas  1994). 

I.  Effectiveness  of  Conversion  Method 

The  conversion  method  developed  executes  in  O(n)  steps,  since  every  amino  acid 
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position  must  be  evaluated.  The  overhead  incurred  by  such  a  conversion  however,  allows 
the  computation  of  the  distances  (and  the  Lennard-Jones  term)  between  all  pairs  of  amino 
acids  in  O (n2)  steps,  since  there  are  n 2  pairs  (length  computations  are  0(1)  in  cartesian 
coordinates).  In  comparison,  computing  the  distance  directly  from  the  protein’s  chemical 
representation  would  incur  an  O (n3)  algorithm,  since  each  length  computation  would 
incur  an  O(n)  penalty.  The  additional  overhead  of  maintaining  a  separate  copy  of  the 
amino  acid  in  cartesian  coordinates  is  more  than  compensated  for  in  decreased  execution 
time. 
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IV.  n -Chain  Hydrocarbons:  A  Sample  Problem 

A.  n-Chain  Hydrocarbon  Problem 

Techniques  aimed  at  solving  the  protein  folding  problem  can  be  tested  on  simpli¬ 
fied  versions  of  the  problem  to  provide  insight  into  the  method’s  particular  characteristics 
and  potential  for  a  solution.  One  such  simplified  version  is  known  as  the  n-chain  hydro¬ 
carbon  problem.  This  problem  attempts  to  find  the  3-dimensional  conformation  of  a  chain 
of  n  hydrocarbons. 

Hydrocarbons  can  be  modeled  in  a  manner  similar  to  proteins:  a  string  of  beads  (a,, 
where  i  =  1 ,...,n),  positioned  by  their  respective  bond  lengths,  bond  angles,  ©,,  and  tor¬ 
sion  angles,  <p(  (see  Chapter  II).  Each  “residue”  in  this  case,  however,  can  be  considered 
identical.  Scientific  knowledge  has  determined  that  all  hydrocarbons  have  “preferred” 
bond  lengths  and  bond  angles.  This  information  can  be  used  in  some  optimization  tech¬ 
niques  to  reduce  the  number  of  problem  parameters.  Note  that  similar  information  has 
been  developed  for  proteins  using  the  Ramachandran  plot  (Lehninger  1970). 

B.  The  Conformation  of  Heptane 

Heptane  is  an  n-chain  hydrocarbon  of  length  seven  (n  =  7).  It  provides  a  simple, 

tractable  problem  for  which  the  conformational  space  can  (almost)  be  exhaustively 
explored.  While  this  structure  reflects  great  simplifications  over  the  protein  folding  prob¬ 
lem,  its  complexity  must  not  be  underestimated. 

Using  the  potential  energy  function  (Chapter  II), 
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Equation  4. 1 


where  r, j  is  the  distance  between  residues,  or  in  this  case  carbons,  separated  by  two  or 

more  residues  along  the  central  chain,  k{  =  310.0  - ^  )>  *e  -  40.0 

Vn  =  1.3  (kcal)  ,  e.  /  =  0.181  (kcal)  ,  o. ,  =  4.0  (Ang) ,  lQ  =  1.526  (Ang)  ,  and 
0O  =  109.47° .  Notice  that  /0  and  0O  represent  the  preferred  bond  angles  and  bond  lengths, 
respectively,  and  n  is  set  to  three,  providing  three  preferred  torsion  angles  at  60°,  180°,  and 
300°,  while  8=0°.  Since  all  bonds  between  hydrocarbons  can  be  considered  identical,  only 
one  value  of  the  ei  t  and  a, j  is  required. 

Even  while  great  simplifications  have  been  made  (modeling  the  molecule  as  a 
chain  of  beads,  assuming  all  residues  to  be  identical,  and  using  a  small  example,  such  as 
heptane),  the  conformation  problem  remains  very  difficult.  Since  torsion  angles  have  three 
preferred  positions  (60°,  180°, 300°),  3”'3,  or  81,  local  minima  are  expected.  (Recall  that 
three  of  the  torsion  angles  are  fixed.)  Tests  on  the  potential  function  support  this  estimate; 
by  starting  from  1000  random  starting  points,  77  local  minima  were  discovered.  However, 
merely  selecting  preferred  angles  does  not  guarantee  that  the  initial  conformer  is  also  a 
local  minimum.  (Notice  in  Table  4.1  that  the  local  minimum  found  by  starting  at  (-60,60,- 


Table  4.1  Selected  Minima  of  Heptane 


m) 

Initial  Position,  <)/0) 

Minimized  Position,  4>* 

-0.339727 

180 

180 

180 

180 

moooooo 

180.000000 

180.000000 

180.000000 

-0,215609 

-60 

60 

-60 

60 

•62.657878 

179.600918 

-179.857154 

179.956953 

-0.215524 

-60 

180 

180 

180 

-62.706316 

179.627701 

180.119720 

180.000000 
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Table  4.1  Selected  Minima  of  Heptane 


u«t>) 

Initial  Position,  0(O) 

Minimized  Position,  0* 

-0.104118 

60 

60 

-60 

60 

62.550498 

180.191468 

-179.798849 

62.592457 

2.097992 

180 

-60 

60 

60 

179.378780 

-79.278332 

79.236470 

62.423103 

4.384072 

180 

-60 

60 

-60 

178.741629 

-93.723491 

69.284953 

-94.158124 

60,60)  results  in  the  same  local  minimum  as  the  point  with  initial  torsion  angles  (- 
60,180,180,180).)  The  global  minimum  for  heptane  corresponds  to  the  point  with  all  tor¬ 
sion  angles  equal  to  180  (see  Table  4.1  and  Figure  4.2).  (For  a  complete  list  of  the  local 


Figure  4.2  Global  Minimum  Conformation,  0*,  for  Heptane 


minima  for  heptane  found  from  the  preferred  torsion  angles,  see  Appendix  A.) 


While  it  may  seem  that  77  local  minima  is  relatively  small,  larger  n-chain  hydro¬ 
carbons  produce  an  exponentially  increasing  number  of  local  minima.  A  molecule  with  as 
few  as  12  beads  has  19,683  possible  preferred  angle  configurations,  with  most  of  these 
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likely  to  produce  a  local  minimum.  It  can  clearly  be  seen  that  finding  the  global  minimum 
for  any  problem  approaching  the  size  of  the  protein  folding  problem  (lengths  of  20  or 
more),  is  intractable  with  exhaustive  methods.  The  next  few  chapters  discuss  methods 
which  attempt  to  find  the  global  solution  without  resorting  to  methods  which  grow  expo¬ 
nentially  intractable  with  problem  size. 
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V.  Convex  Global  Underestimator 

A.  A  Global  Underestimator 

One  possibility  for  aiding  the  search  for  the  global  minimum  of  the  molecule’s 
potential  energy  function  is  to  use  a  global  underestimator  to  localize  the  search  in  the 
region  of  the  global  minimum.  This  new  method  is  designed  to  fit  all  known  local  minima 
with  a  convex  function  which  underestimates  all  of  them,  but  which  differs  from  them  by 
the  minimum  possible  amount  in  the  discrete  Li  norm.  The  minimum  of  this  underestima¬ 
tor  is  used  to  predict  the  global  minimum  for  the  function,  allowing  a  localized  conformer 
search  to  be  performed  based  on  the  predicted  minimum.  A  new  set  of  conformers  gener¬ 
ated  by  the  localized  search  can  then  serve  as  a  basis  for  another  quadratic  underestima¬ 
tion.  After  several  repetitions,  the  global  minimum  can  be  found  with  reasonable 
assurance. 

The  use  of  an  underestimating  function  allows  the  translation  of  a  very  complex 
function  into  a  simple  underestimator.  If  the  underestimator  is  well-suited  to  the  problem 
(provides  accurate  predictions  for  the  global  minimum),  immense  savings  in  time  can  be 
accomplished.  The  presence  of  quadratic  terms  in  the  length  and  bond  angle  portions  of 
the  molecular  energy  function  (See  Equation  2.1.)  indicate  the  use  of  a  convex  quadratic 
may  provide  a  suitable  approximation  for  this  problem. 

B.  Theory 

This  method  is  presented  in  terms  of  a  differentiable  function  F(<f>),  where  4>  e  Rx 
(t  represents  the  number  of  torsion  angles),  and  where  F(<}>)  has  many  local  minima.  Since 
the  bond  lengths  and  bond  angle  terms  in  the  potential  function  carry  severe  quadratic 
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penalties,  they  can  be  assumed  to  be  fixed  for  the  purposes  of  the  global  underestimator. 
Therefore,  F(<J>)  is  a  vector  of  torsion  angles  of  length  x,  where  x  =  n  -  3  and  n  is  the  num¬ 
ber  of  residues.  (Recall  that  the  first  three  are  fixed  by  definition.)  To  begin  the  iterative 
process,  a  set  of  k  >  2x+l  distinct  local  minima  are  computed.  This  can  be  done  with  rela¬ 
tive  ease  by  using  an  efficient  unconstrained  minimizer,  starting  with  a  large  set  of  points 
chosen  at  random  in  the  hypercube  H<f>,  which  is  assumed  to  enclose  the  torsion  angle 
space. 

Assuming  that  k  >  2x+l  local  minima  0^,  for  j=l,...,k,  have  been  computed,  a  con¬ 
vex  quadratic  underestimator  function  'F(<j>)  is  now  fitted  to  these  local  minima  so  that  it 
underestimates  all  the  local  minima  and  normally  interpolates  F(0^)  at  2x+l  points.  This 
is  accomplished  by  determining  the  coefficients  in  the  function  4i(0)  so  that: 

6 j  =  F  (0  U) )  -  (0  O') )  >0  Equation  5. 1 

for  j=l,...,k,  and  where  { 5^ is  minimized.  That  is,  the  difference  between  F(0)  and 

VF(0)  is  minimized  in  the  discrete  Li  norm  over  the  set  of  k  local  minima  0^,  j=l . k.  The 

underestimating  function  4^(0)  is  given  by: 

i  ,  . 

¥(0)  =  co  +  X  (  CI4I  +  \dM)  •  Equation  5.2 

i  =  1 

Note  that  Cj  and  dj  appear  linearly  in  the  constraints  (Equation  5.1)  for  each  local 
minimum  0^.  Convexity  of  this  quadratic  function  is  guaranteed  by  requiring  that  d,  >  0 
for  i-l,...,x.  Other  linear  combinations  of  convex  functions  could  also  be  used,  but  this 


quadratic  function  is  the  simplest. 
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C.  Minimizing  the  Underestimating  Function 

The  unknown  coefficients  ci(  i=0,...,T,  and  dr  i=l,...,x,  can  be  determined  by  a  lin¬ 
ear  program  which  may  be  considered  to  be  in  the  dual  form.  For  reasons  of  efficiency,  the 
equivalent  primal  of  this  problem  is  actually  solved,  as  described  below.  The  solution  to 
this  primal  linear  program  provides  an  optimal  dual  vector,  which  immediately  gives  the 
underestimating  function  coefficients  q  and  dj.  Since  the  convex  quadratic  function  ¥(<(>) 
gives  a  global  approximation  to  the  local  minima  of  F(<J>),  then  its  easily  computed  global 
minimum  function  value  is  a  good  candidate  for  an  approximation  to  the  global  min¬ 
imum  of  the  correct  energy  function  F(<j>). 

An  efficient  linear  programming  formulation  and  solution  satisfying  Equation  5.1 
and  Equation  5.2  will  now  be  summarized.  Let  =  F(<{)^),  for  j=l,...,k,  and  let  /  e  Rk  be 
the  vector  with  elements  f®.  Also  let  to®  e  RT  be  the  vector  with  elements  1  /2(<h(j)I)2, 
i=l,...,T,  and  ek  €  Rk  be  the  vector  of  ones.  Now  define  the  following  two  matrices  <f>  e 
rCc+OxIc  and  Q  €  R^: 


Equation  5.3 


Finally,  let  c  e  RT+1,  d  €  RT,  and  8  e  Rk  be  the  vectors  with  elements  q,  dp  and  8,,  respec¬ 
tively.  Then  (Equation  5.1)  and  (Equation  5.2)  can  be  restated  as  the  linear  program  (with 
variables  c,  d,  and  8): 


minimize  eTZ 
c,  d,  8 
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<Dr  QT  0 
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r/> 

j 

subject  to 

-<&T  -QT  -Ik 
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-f 
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o 

i 

o 

8 

_0_ 

Equation  5.4 


It  is  easy  to  see  that  c  =  d  =  0  and  8  =  /  >  0  is  always  feasible  for  (Equation  5.4), 
and  since  efc78  >  0,  the  linear  program  always  has  an  optimal  solution.  Furthermore,  since 
the  matrix  in  (Equation  5.4)  has  more  rows  than  columns  (2k+x  rows  and  k+2x+l  col¬ 
umns,  where  k  >  x+1),  it  is  computationally  more  efficient  to  consider  it  as  a  dual  problem, 
and  to  solve  the  equivalent  primal.  After  some  simple  transformations,  this  primal  prob¬ 
lem  reduces  to: 


minimize 


subject  to 


0>  0 
Q  -ix\ 


fTy\-fTek 


y\ 

®ek 

?ek 

ylfy2^0 


Equation  5.5 


which  has  only  2x+l  rows  and  k+x  >  3x+l  columns,  and  the  obvious  initial  feasible  solu¬ 
tion  yj  =  e^  and  y2  =  0.  Its  optimal  solution  gives  the  values  of  c,  d,  and  8  via  the  dual  vec¬ 
tors,  and  also  determines  which  values  of  are  interpolated  by  the  potential  function 
'FOJ)).  That  is,  the  basic  columns  in  the  optimal  solution  to  (Equation  5.5)  correspond  to 
the  conformations  for  which  =  4/(4><J-)). 

Note  that  once  an  optimal  solution  to  (Equation  5.5)  has  been  obtained,  the  addi¬ 
tion  of  new  local  minima  is  very  easy.  It  is  done  by  simply  adding  new  columns  to  <I>  and 
Q  and  therefore  to  the  constraint  matrix  in  (Equation  5.5).  The  number  of  primal  rows 
remains  fixed  at  2x+l,  independent  of  the  number  k  of  local  minima. 
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D.  Iteration  -  Approaching  the  Global  Minimum 

The  convex  quadratic  underestimating  function  4/(<j>)  determined  by  the  values  c  e 

Rx+1  and  d  €  RT  now  provides  a  global  approximation  to  the  local  minima  of  F(<b),  and  its 
easily  computed  global  minimum  point  <})min  is  given  by  (<bmin)i  =  i=l»—J,  with  cor¬ 

responding  function  value  given  by  vPm]n  =  cQ  -  £}=  ]  c}/  (2 d.)  .  The  value  'Fmin 
is  a  good  candidate  for  an  approximation  to  the  global  minimum  of  the  correct  energy 
function  F(<|)),  and  therefore  0min  can  be  used  as  an  initial  starting  point  around  which 
additional  configurations  (i.e.  local  minima)  should  be  generated.  These  local  minima  are 
added  to  the  constraint  matrix  in  (Equation  5.5)  and  the  process  is  repeated.  Before  each 
iteration  of  this  process,  it  may  be  advantageous  to  reduce  the  volume  of  the  hypercube 
H<fi  over  which  the  new  configurations  are  produced  so  that  a  tighter  fit  of  ^((J))  to  the  local 
minima  “near”  may  be  produced. 

The  rate  and  method  by  which  the  hypercube  size  is  decreased  and  the  number  of 
additional  local  minima  computed  at  each  iteration  will  be  determined  by  computational 
testing.  The  selection  of  the  hypercube  size  can  be  done  computationally  or  with  the  aid  of 
an  expert.  But  clearly  the  method  depends  most  heavily  on  computing  local  minima 
quickly  and  on  solving  the  resulting  linear  program  efficiently  to  determine  the  approxi¬ 
mating  function  4/(<j>)  over  the  current  hypercube. 

If  Ec  is  a  cutoff  energy,  then  one  means  for  decreasing  the  size  of  the  hypercube 
H<>  at  any  step  is  to  let  H<j)  =  {<(>:  H^)  <  Ec).  To  get  the  bounds  of  H<j>,  consider  H/(4>)  <  Ec 
where  VF(<}>)  satisfies  (Equation  5.2).  Then  limiting  <)>;  requires  that: 

l(cA  +  ^?)s*e-cO- 

i  =  1 


Equation  5.6 


33 


As  before,  the  minimum  value  of  ¥(<)>)  is  attained  when  0,  =  -c/dj,  i=l,...,T. 
Assigning  this  minimum  value  to  each  0P  except  0k,  then  results  in: 

c*0*  +  \dfil  <  Ec  -  c0  +  \  £  cf/di  •  .  Equation  5.7 

Z  Zi*k 

The  lower  and  upper  bounds  on  0^,  k=l,...,t,  are  given  by  the  roots  of  the  quadratic 
equation 

ck$k  +  \dk$k  ~  P/t  •  Equation  5.8 

Hence,  these  bounds  can  be  used  to  define  the  new  hypercube  H0  in  which  to  gen¬ 
erate  new  configurations. 

E.  Local  Minimization 

The  final  predicted  point  can  be  used  as  the  initial  point  for  a  local  minimization 
procedure.  Since  the  point  predicted  by  the  underestimator  is  expected  to  be  very  close  to 
the  global  minimum,  the  conformation  found  at  the  nearest  local  minimum  should  corre¬ 
spond  to  the  global  minimum. 

An  accurate  and  efficient  local  minimization  procedure  is  of  great  importance  to 
the  establishment  of  an  accurate  underestimator  as  well.  The  effectiveness  of  the  global 
underestimating  procedure  relies  heavily  upon  the  ability  to  produce  very  many  local  min¬ 
ima.  Since  this  minimization  is  performed  a  large  number  of  times  and  each  minimum 
must  reflect  an  accurate  conformer,  the  development  an  efficient  algorithm  is  extremely 


important. 
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VI.  Finding  Local  Minima 

A.  Twice  Differentiable  Unconstrained  Optimization 

Since  the  potential  energy  function  is  smooth  (twice-continuously  differentiable), 

any  method  which  computes  an  n-dimensional  unconstrained  minimum  can  be  used.  The 
methods  to  be  examined  satisfy  the  basic  form  (where  tyk  is  the  current  estimate  of  <}>*,  the 
global  minimizer  conformation): 

1 .  Obtain  some  initial  conformation  cp0  and  set  k  «-  0 . 

2.  [Test  for  convergence.]  If  the  conditions  for  convergence  are  satisfied,  the 
algorithm  terminates  with  tyk  as  the  solution. 

3.  [Compute  a  search  direction.]  Compute  a  non-zero  n -vector  pk,  the  direction 
of  search. 

4.  [Compute  a  step  length.]  Compute  a  positive  scalar  ak,  the  step  length,  which 

satisfies  the  descent  condition,  F  +  a kpk)  <  F(§k) . 

5.  [Update  the  estimate  of  the  minimum.]  Set  §k  +  ,  <-  +  a kpk ,  k  <—  it  +  1 ,  and 

go  back  to  step  2. 

To  satisfy  the  descent  condition  (that  for  a  positive  scalar,  F(< \>k  +  o.kpk)  <  F(< pk)  ), 
pk  and  ak  must  have  certain  properties.  One  way  of  guaranteeing  that  F  can  be  reduced  at 
the  kth  iteration  is  to  require  that  pk  satisfies 

T 

gkPk  <  0  Equation  6. 1 

where  gk  =  VF(< fyk)  is  the  gradient  of  the  function,  F,  at  point  <()*.  When  this  condition  is 
satisfied,  pk  is  known  as  a  descent  direction.  The  next  three  sections  will  describe  three 
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related,  but  different,  methods  for  computing  the  local  minima  of  F(<j>). 


B.  Steepest  Descent 

Perhaps  the  greatest  impact  on  the  effectiveness  of  a  local  minimization  algorithm 
is  the  selection  of  a  “good”  search  direction,  pk.  Even  in  as  few  as  two  dimensions  there 
are  an  infinite  number  of  ways  to  compute  the  descent  direction.  The  importance  of  this 
decision  led  to  the  naming  of  local  minimizers  according  to  the  procedure  by  which  it 
computes  the  search  direction. 

Since  a  linear  approximation  to  the  function  F(<J>)  can  be  presented  as  a  Taylor 
series  expansion  about  xk  (where  gk  denotes  g  (§k)  =  VF(< j>p  ), 

T 

F(tyk  +  p)  =F(tyk)  +  gkp  Equation  6.2 


T 

it  would  appear  that  a  large  reduction  in  F  would  occur  when  gkp  is  large  and  negative. 
However,  limits  must  be  placed  on  p.  Otherwise,  for  any  p  such  that  gkp  <  0  ,p  could  sim¬ 
ply  be  chosen  as  an  arbitrarily  large  multiple  of  p .  The  aim  is  to  chose  p,  so  that  among  all 

J 

normalized  vectors,  gkp  is  a  minimum.Therefore,  given  some  norm,  HI  ,pk  is  the  solution 
of  the  minimization  problem: 


minimize 

P£3!" 
bl  =  i 


[SkP 


)■ 


Equation  6.3 


If  the  two-norm  is  used  for  the  calculation  of  ||/?|| ,  the  solution  to  Equation  6.3  is  just  the 
negative  gradient, 

Pk  =  ~g  (4>*)  •  Equation  6.4 


This  solution  is  known  as  the  direction  of  steepest  descent.  Clearly,  as  long  as  the  gradient 
does  not  equal  zero,  a  descent  direction  is  guaranteed. 
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Unfortunately,  while  a  steepest  descent  algorithm  is  relatively  simple  and  is  guar¬ 
anteed  to  converge  (reach  a  minimum),  it  is  not  very  efficient.  When  faced  with  a  function 
of  quadratic  or  greater  complexity,  the  steepest  descent  algorithm  may  require  several 
hundred  iterations  to  make  very  little  progress  towards  the  solution.  The  function  F(<j>)  pre¬ 
sents  just  such  a  problem.  Hence  more  efficient  methods  must  be  explored  (Gill  et  al. 
1981). 

C.  Newton’s  Method 

The  primary  weakness  of  the  steepest  descent  method  is  its  attempt  to  represent  a 
complex  function  using  a  linear  model  involving  only  the  gradient.  The  use  of  a  quadratic 
model  retains  much  of  the  simplicity  found  in  the  steepest  descent  procedure,  yet  can  yield 
greater  success  and  efficiency  in  practice. 

This  model,  based  on  the  first  and  second  derivatives  of  the  function  F,  is  known  as 
Newton’s  Method.  The  function  F  is  approximated  using  the  first  three  terms  of  the  Tay- 
lor-series  expansion: 

F  ( xk  +p)  ~  Fk  +  g[p  +  \pTGkp  Equation  6.5 


2 

where  G*  denotes  the  Hessian,  V  F(<J>)  at  the  point  tyk.  It  is  helpful  to  formulate  the  algo¬ 
rithm  in  terms  of  p  (the  step  to  the  minimum),  rather  than  the  predicted  minimum  itself.  A 
minimum  of  the  right  hand  side  of  Equation  6.5  will  be  achieved  at  the  minimum  of  the 
quadratic: 


O(p) 


T  1  T 

=  8kP  +  iP  GkP 


Equation  6.6 


The  function  O  has  a  stationary  point  only  if  there  exists  a  point  where  its  gradient  vector 
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vanishes  (becomes  zero).  Therefore,  for  any  minimum,  the  gradient  of  <l>  must  be  equal  to 
zero;  that  is,  VO(p)  =  Gkp  +  gk  =  0.  A  minimization  algorithm  in  which  pk  is  defined 
by: 

pk  =  ~Gkgk  Equation  6.7 

is  termed  Newton’s  method,  and  the  direction,  pk ,  a  Newton  direction. 

If  Gk  is  positive  definite  and  F(0)  is  itself  quadratic,  only  one  iteration  is  required 
to  reach  the  minimum  of  the  model  function  from  any  starting  point.  Therefore,  good  con¬ 
vergence  from  the  Newton’s  method  is  expected  when  the  quadratic  model  is  accurate.  For 
a  general  function  (F)  Newton’s  method  converges  quadratically  to  0*  if  the  starting  point 
(<(>0)  is  sufficiently  close  to  0*,  the  Hessian  matrix  is  positive  definite  at  0*,  and  the  step 
lengths  ( ak )  converge  to  unity. 

The  local  convergence  properties  of  Newton’s  method  make  it  an  exceptionally 
attractive  algorithm  for  unconstrained  optimization.  The  availability  of  the  second  deriva¬ 
tive  also  allows  the  verification  that  sufficient  conditions  for  optimality  exist.  (Gill  et  al. 
1981)  While  Newton’s  method  can  be  an  extremely  effective  optimization  method,  it 
requires  the  often  very  difficult  computation  of  second  derivatives. 

D.  Quasi-Newton  Methods 

Indeed,  F(0)  presents  a  very  difficult  case  since  not  only  are  the  second  derivatives 
difficult  to  compute,  so  are  the  first!  Ideally,  advantage  of  the  curvature  information  (the 
Hessian  matrix)  could  be  taken  without  having  to  compute  it  directly.  Just  such  an  attempt 
is  made  in  quasi-Newton  methods.  These  algorithms  are  based  on  the  idea  of  “building 
up”  curvature  information  as  the  iterations  of  a  descent  method  proceed.  Thus,  an  approx- 
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imation  to  the  curvature  of  a  non-linear  function  can  be  computed  without  explicitly  form¬ 
ing  the  Hessian  matrix. 

Quasi-Newton  methods  are  very  similar  to  Newton’s  method,  except  that  G^'is 
approximated  by  a  symmetric,  positive  definite  matrix,  ,  which  is  corrected  or  updated 
from  iteration  to  iteration.  The  basic  structure  for  the  kth  iteration  is: 

1 .  Set  the  search  direction,  = 

2.  Compute  a*  using  line  search  along  giving  =  <\>k  +  aksk . 

3.  Update  H*  giving  H*+1. 

The  initial  matrix,  Hj  can  be  any  positive  definite  matrix,  although  in  the  absence  of  better 
estimates,  Hj=I  is  often  made. 

The  update  of  the  Hessian  matrix  is  generally  accomplished  through  the  use  of  an 
updating  formula.  Many  formulas  have  been  established,  each  of  which  attempts  to  gener¬ 
ate  a  Hessian,  H^+1,  from  the  previous  Hessian  H*  and  preserve  positive  definiteness. 
Through  the  use  of  second  derivative  information  gained  on  the  kth  iteration,  repeated 
updating  can  change  the  arbitrary  matrix,  Hj,  into  a  close  approximation,  H*,  to  . 

One  possibility  for  updating  the  Hessian  was  developed  by  Broyden,  Fletcher, 
Goldfarb,  and  Shanno,  and  is  known  as  the  BFGS  formula  (Fletcher  1987).  If  the  differ¬ 
ences,: 

6 k  =  a ksk  =  tyk+l-tyk  Equation  6.8 

yk  =  8k  +  ,  -  gk  Equation  6.9 

are  defined,  then  the  BFGS  formula  can  be  stated  as: 
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H 
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(  T  \ 

55^_ 

'  8yrH  +  Hy8r 

K  bTy  y 

8ry 

k  8ry  > 

Equation  6.10 


Note  that  the  subscript,  k,  has  been  suppressed  on  the  right  hand  side.  While  many  other 
formulas  have  been  developed,  none  have  shown  any  significant  improvement  over  the 
BFGS  method. 

Quasi-Newton  methods  provide  several  distinct  advantages  over  a  Newton’s 
method  optimization: 

1 .  Only  first  derivatives  are  required. 

2.  Positive  definiteness  of  H*  will  guarantee  the  descent  property. 

3.  Only  0(n2)  multiplications  per  iteration  are  required.  (0(n3)  are  required  for 
Newton’s  method.) 

E.  Finite  Difference  Computations 

Unfortunately,  computation  of  even  the  first  derivative  can  be  a  source  of  great  dif¬ 
ficulty,  as  is  the  case  with  the  potential  function  given  by  due  to  the  mixture  of  cartesian 
and  chemical  representations.  However,  the  ability  to  expand  smooth  functions  in  Taylor 
series  allows  the  values  of  derivatives  to  be  approximated.  This  is,  if  hj  is  the  finite  differ¬ 
ence  interval  (the  spacing  over  which  the  approximation  is  to  be  made)  associated  with  the 
j- th  component  of  <j>  and  e;  is  the  vector  with  unity  in  the  j- th  position  and  zeros  elsewhere, 
then: 

FiQ  +  hjej)  =  F  (<j> )  +  hjg Tej  +  ^ hjejG  (<£>  +  ?, hjej)  e-  Equation  6. 1 1 

T  T 

where  0  <  p  1  <  1 .  Since  the  terms  g  e.  and  e.  G  (0  +  p , h.e.)  are  the  y'-th  element  of  g, 
and  the  j-t h  diagonal  element  of  G  (()>  +  p,/yr)  ,  respectively,  then: 
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gj  =  -j-  (F  (<}>  +  hjej)  -  F  (<j>) )  +  0  ( hj )  .  Equation  6. 12 

Equation  6.12  describes  the  forward-difference  interval  expression  in  the  direction  ey  The 
error  developed  in  approximating  gj  by  the  expression  to  the  left  of  Equation  6.12  is  due  to 

l  2  r 

neglecting  the  term,  -hj  G (0  +  Pj/rtr)  e}.  This  inaccuracy  is  known  as  the  truncation 
error ,  since  it  arises  from  the  truncation  of  the  Taylor  series.  In  general,  p!  is  unknown  and 
therefore  it  is  possible  only  to  compute  or  estimate  an  upper  bound  on  the  truncation  error. 
Similarly,  a  backward-difference  expression  can  be  developed  by  expanding  the  Taylor 
series  about  the  point  (6  -  ft.e.)  in  order  to  get: 

gj  =  ^  (F  (40  -  F  (4>  -  hjCj) )  +  O  (hj)  .  Equation  6. 13 

If  the  Taylor  series  expansions  of  F(4>)  using  both  the  forward  and  backward  difference 
equations  are  carried  one  term  further,  the  result  is: 

F  (4)  +  hjCj)  =  F  (40  +  hjgjej  +  ^h^e^Gjej  +  0^  H*  j  Equation  6. 14 

and: 

FW-hjej)  =  F(40  ~hJgJej  +  ^ hj  ej G]ej  +  0^  hj  j .  Equation  6. 1 5 

Subtracting  and  dividing  by  hj,  one  obtains: 

£,  (40  =  +  V/)  ~FW ~jjej) ]  +  h] )  Equation  6. 1 6 

since  the  terms  containing  G,  cancel.  The  right  hand  side  of  this  equation  is  known  as  a 
central-difference  approximation  to  £,(4>).  This  formulation  provides  a  truncation  error  of 
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the  second  order,  involving  the  third  derivative  of  F  in  the  region  [ty-hjej,  ty+h^j].  How¬ 
ever,  two  function  evaluations  (in  addition  to  F(<j>))  are  required  to  approximate  gj(§),  as 
compared  to  one  evaluation  for  Equation  6.12  and  Equation  6.13  (the  forward  and  back- 
ward-difference  approximations). 

Through  the  use  of  the  central  difference  approximation  (Equation  6.16),  the  gra¬ 
dient  of  the  function  can  be  calculated  without  explicit  knowledge  of  the  derivative.  The 
error  produced  in  the  finite  difference  approximation  should  be  small  as  long  as  hj  is  cho¬ 
sen  to  be  suitably  small  (Gill  et  al.  1981). 

F.  Step  Length  (Line  Search  Calculations) 

Once  the  search  direction,  s,  has  been  established,  the  distance  to  travel  along  this 

path,  a,  must  be  found.  Any  line  search  method  will  suffice  for  this  calculation.  Since  the 
potential  function  is  of  unknown  complexity,  a  simple  algorithm  is  chosen: 

1.  Let/0  =  F«j». 

2.  Set  oco  to  a  very  small  initial  value  (a  «-  e  >  0 )  and  i  <-  0 . 

3.  Compute  a  new  function  value,  f.  j  =  F(<t>  +  a.s) . 

4.  If  /l  + 1  <  ft ,  set  a.  +  j  =  2a(  ,/<-/+  1 ,  and  return  to  step  3. 


5. 


Otherwise,  set  a  =  a  and  stop. 
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VII.  Underestimating  Function-  Trials  and  Results 

A.  Global  Minimum  Conformation  for  Heptane 

A  solution  for  the  7  bead  hydrocarbon  heptane  was  obtained  in  68  seconds  of  CPU 

time  on  a  SUN  Sparcstation  2.  While  the  improvements  in  computation  time  are  small  for 
this  problem  (Enumeration  required  only  79  seconds  to  test  all  local  minima.),  many  of  the 
characteristics  of  the  method  can  be  observed.  Table  7.1  summarizes  the  minimization  of 
the  potential  function  for  heptane.  (See  Appendix  C  for  a  complete  listing  of  the  program 
output.)  Notice  that  the  first  iteration  reports  the  current  “best”  (corresponding  to  the  low¬ 
est  potential  energy  value)  conformation  to  be  approximately  the  set  of  torsion  angles, 
(180\180\300\180°).  Recall  that  the  bond  lengths  and  bond  angles  are  considered  fixed 
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for  the  purposes  of  the  global  underestimator.  However,  the  global  underestimator  predicts 
a  conformation  close  to  (180°,  180°,  180°,  180°).  After  the  second  iteration,  concentrated 
around  the  predicted  point,  the  new  results  report  that  a  better  function  value  can  be  real¬ 
ized  at  this  conformation.  The  underestimator  and  the  current  global  minimum  now  reside 
at  almost  identical  points.  The  third  iteration  confirms  (180°,  180°,  180°,  180°)  as  the  glo¬ 
bal  minimum.  This  result  has  been  shown  through  enumeration  to  be  the  correct  global 
minimum  conformation  for  the  potential  energy  function. 

B.  Prediction  for  Hydrocarbons  of  Arbitrary  Length 

As  problem  size  increases,  the  number  of  local  minima  increases  exponentially 

(0(3n)),  as  can  be  seen  in  Figure  7.2.  With  exponential  increases  in  the  number  of  local 


Number  of  Residues 

Figure  7.2  Number  of  Residues  vs.  Number  of  Local  Minima 

minima,  exhaustive  searches  for  the  global  minimum  quickly  become  computationally 
intractable.  By  comparison,  the  underestimator  method  produces  a  final  prediction  in 
O (n4)  time. 

The  main  contributors  to  the  underestimator’s  O (n4)  rate  are  the  function  evalua- 
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tions  and  conversions.  Every  function  evaluation  incurs  an  0(n2)  penalty  due  to  the  con¬ 
versions  and  distance  calculations  required  by  the  Lennard-Jones  term.  The  non-linear 
program  which  computes  the  local  minima  averages  O(n)  iterations  to  reach  a  local  mini¬ 
mum  (Figure  7.3).  The  O(n)  iterations  performed  in  the  local  minimization  procedure. 


Figure  7.3  Number  of  Iterations  Required  for  Local 
Minimization 

combined  with  the  O (n2)  time  for  function  evaluation,  indicate  that  the  total  time  required 
per  local  minimization  increases  at  a  rate  no  greater  than  0(n  )  (See  Figure  7.4.).  In  fact, 
the  time  required  to  perform  a  single  local  minimization  as  a  function  of  the  number  of 
beads,  n,  can  be  approximated  (using  a  least  squares  method)  by  the  function 
f (n)  =  0.003 n  -  0.68 n  +  5.46n  -  14.28  (as  seen  by  the  curve  in  Figure  7.4). 

The  amount  of  time  for  exhaustive  techniques  is  driven  by  the  number  of  local 
minimizations  required  to  test  each  case  (531,441  minimizations  in  the  15  bead  problem), 
and  hence,  to  test  every  possible  local  minimum  conformation  for  an  n-chain  hydrocarbon 
of  length  15  would  take  over  185  days  on  a  Sun  workstation!  In  contrast  to  enumeration, 
the  underestimator  requires  a  much  less  drastic  increase  in  the  number  of  local  minima 
calculations.  Since  the  number  of  iterations  of  the  underestimator  needed  to  find  a  global 
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Figure  7.4  Time  for  Local  Minimization  per  Number  of 

Residues 


minimum  remains  close  to  constant  (0(1)),  and  each  step  requires  O(n)  local  minima  to 
form  the  underestimator,  a  linear  increase  in  the  number  of  minimizations  is  expected. 

This  effect  can  be  seen  in  Figure  7.5. 
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Figure  7.5  Number  of  Local  Minima  for  Complete  Global 
Minimum  Prediction 

The  effect  of  increases  in  problem  size  on  total  time  for  global  minimum  predic¬ 
tion  can  be  seen  in  Figure  7.6.  Note  that  this  graph  depicts  CPU  time  and  does  not  include 
the  time  required  to  enter  values  between  iterations.  The  curve  seen  in  Figure  7.6  repre- 
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Figure  7.6  Total  CPU  Time  for  Global  Minimum  Prediction 

4  3  2 

sents  the  least  squares  fit,  0.01 «  -  0.2 9n  +  3.85  n  -  22.26 n  +  44.13 .  Thus  the  final  predic¬ 
tion  can  be  provided  in  no  worse  than  O (n4)  time.  While  the  prediction  times  for  the 
underestimator  are  clearly  significantly  better  than  those  of  an  exhaustive  search,  the 
O (n4)  rate  of  increase  in  time  as  problem  size  grows  makes  it  unlikely  that  molecules  sub¬ 
stantially  larger  than  22  beads  will  be  solved.  In  fact,  the  underestimating  method  enabled 
predictions  to  be  made  for  hydrocarbons  of  length  4  through  15  and  22.  (See  Appendices 

B,  C,  D,  and  E  for  selected  underestimator  runs.) 

The  22  bead  problem  took  over  16  hours  for  the  final  prediction  (Figure  7.7).  This 
prediction  was  made  with  4  iterations  of  the  global  underestimator  (See  Appendix  E  for 
complete  output  from  the  22  bead  run.),  and  143  seconds  were  required  for  each  local 
minimization.  Further  increases  in  problem  size  will  face  even  larger  increases  in  execu¬ 
tion  time. 

C.  Applicability  to  the  Protein  Folding  Problem 

While  the  molecular  conformation  problem  is  a  greatly  simplified  version  of  the 
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Figure  7.7  Global  Energy  Minimum  for  a  22 
bead  hydrocarbon 


protein  folding  problem,  the  basic  “unified  atom”  model  used  is  still  applicable  to  both. 
The  only  substantial  change  is  that  the  potential  energy  function  should  be  appropriate  to 
the  protein  structure.  The  new  potential  energy  function  would  be  required  to  model  inter¬ 
actions  between  all  combinations  of  the  20  amino  acids,  as  well  as  take  into  account  any 
additional  potential  energy  contributions  caused  by  amino  acid  interactions.  Notice,  how¬ 
ever,  that  this  model  completely  ignores  any  intra-peptide  interactions;  that  is,  the  model 
used  does  not  consider  any  secondary  structures  such  as  alpha  helixes  and  beta  strands.  A 
much  more  complicated  model  would  be  required  to  capture  these  interactions.  Neverthe¬ 
less,  the  formulation  of  the  convex  quadratic  global  underestimator  and  the  local  minimi¬ 
zation  procedures  based  on  the  model,  presented  in  Chapter  II,  will  remain  applicable  to 
the  problem. 


VIII.  Quadratic  Assignment  Formulation 


A.  Lattice  Restrictions 

Attempts  to  solve  the  protein  folding  problem  in  continuous  space  confront  the 
daunting  task  of  finding  the  optimal  solution  from  among  an  infinite  number  of  conforma¬ 
tions.  Lattice  methods  attempt  to  avoid  this  by  restricting  each  amino  acid  to  a  vertex,  or 
site,  on  a  fixed  lattice.  By  limiting  the  positions  of  the  residues,  the  potential  energy  mini¬ 
mization  can  be  formulated  in  a  manner  which  facilitates  the  optimization.  This  section 
describes  a  method  which  involves  a  quadratic  assignment  formulation  for  the  protein 
folding  problem  using  the  lattice  model.  The  quadratic  assignment  problem  can  then  be 
transformed  into  a  continuous  concave  quadratic  global  minimization  problem.  The  global 
solution  to  this  concave  minimization  problem  can  then  be  used  as  a  starting  point  for  a 
“relaxed”  continuous  minimization  problem.  The  result  this  minimization  should  provide 
a  global,  or  near  global,  minimum  of  the  potential  energy  function,  and  hence  a  prediction 
of  the  native,  or  folded  state,  of  the  molecule  (Phillips  and  Rosen  1994). 

B.  Quadratic  Assignment  Formulation 

To  formulate  a  discrete  approximation  to  the  molecular  conformation  problem,  the 

original  continuous  problem  in  3-dimensional  space  is  approximated  by  a  suitable  3- 
dimensional  lattice  with  N  sites,  N>n  (where  n  is  the  number  of  amino  acids).  If  Sj  repre¬ 
sents  lattice  site  j,  for  j=  1,...,jV,  then  a  total  of  nN  zero-one  variables  x,j  are  sufficient  to 
completely  determine  the  assignment  of  the  beads  a,  to  the  lattice  sites  sj.  Two  constraints 
are  established: 

1 .  Each  bead  must  occupy  exactly  one  lattice  site. 
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i  =  1,  ...,n 


Equation  8. 1 


2.  At  most  one  bead  occupies  each  lattice  site,  sj. 

n 

XX«fS1'  '  *  i,—,N 

j  -  1  Equation  8.2 

The  objective  function  can  be  written  in  two  portions,  a  linear  and  quadratic  term.  The  lin¬ 
ear  term: 

n  N 

X  X  dijxi]  Equation  8.3 

«=U=1 

represents  the  direct  contribution,  dy,  to  the  potential  energy  when  the  bead  a,  is  assigned 
to  lattice  site  sj.  Direct  contributions  consist  of  the  interaction  of  the  residue,  a;,  and  the 
environment  (hydrophobic  beads  might  prefer  to  be  in  the  interior  of  the  folded  chain  and 
hence  in  the  center  of  the  lattice  structure).  The  quadratic  term: 

j  n  N  n  N 

9XI  X  X  Pijkixijxkl  Equation  8.4 

» =  1  j  =  1  *  =  1  /  =  1 

represents  the  pairwise  contribution  to  the  total  potential  energy  when  the  bead  %  is 
assigned  to  lattice  site  Sj  and  bead  a*  is  assigned  to  lattice  site  s/.  One  possibility  for  this 
penalty  is,  as  in  the  continuous  case,  the  Lennard-Jones  pairwise  potential.  Using  Len- 
nard-Jones,  the  term  pijki  depends  only  on  a,,  a*,  and  I  Is,  -  s/ll.  Therefore,  it  is  convenient  to 
write  the  total  potential  energy  function  in  the  form: 

E  (jc)  =  c  x  +  ^  x  Qx 


Equation  8.5 
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where  x  e  RnN  denotes  the  zero-one  vector  with  elements  xip  c  e  RnN  denotes  the  vector 
with  elements  dip  and  Qe  R(n/V)  x  is  a  real  symmetric  (typically  indefinite)  matrix 
with  elements  p^.  Thus,  the  3-dimensional  minimum  energy  conformation  of  the  lattice 
restricted  model  is  given  by  the  solution  to  the  quadratic  assignment  problem : 

min  E  ( x ) 


Equation  8.6 


where  P  consists  of  the  constraints  (Equation  8.1)  and  (Equation  8.2),  and  the  integer 
restrictions,  x^  e  (Q 1}  for  i=  1 and  j-  1 

The  usefulness  of  this  approach  is  very  heavily  dependent  on  the  selection  of  the 
lattice  structure.  Obviously,  determining  the  perfect  lattice  is  as  difficult  as  the  protein 
folding  problem  itself.  Great  care  must  be  taken  to  insure  the  lattice  is  fine  enough  to 
closely  approximate  possible  solutions  in  continuous  space,  and  coarse  enough  so  that  the 
problem  does  not  become  computationally  intractable.  The  lattice  structure  can  take  many 
forms:  spherical,  rectangular,  random,  or  specially  developed  (Hinds  and  Levitt  1992). 
That  the  lattice  might  not  allow  for  the  global  minimum,  or  any  stable  conformation  at  all, 
is  not  important.  The  continuous  minimization  stage  will  relax  the  “non-stable”  conformer 
found  in  the  lattice  restricted  phase  to  a  true  energy  minimum. 


C.  Concave  Quadratic  Global  Minimization  Formulation 

The  discrete  quadratic  assignment  problem  (Equation  8.6)  can  be  shown  to  be 

equivalent  to  the  continuous  minimization  of  a  strictly  concave  quadratic  function  over  a 

polytope.  If  Xmax  is  the  maximum  eigenvalue  of  the  matrix  Q  =  ,  and 

2 

p  =  1  +Xmajf,  then  since  =  xi}  (recall  that  x(j  e  {Q 1}  ),  the  energy  function  can  be 


rewritten  as: 
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E'  (x)  =  c,T x  +  ^xTQ'x 


Equation  8.7 


where  c  =  c  +  [^je  and  Q'  =  Q-\il  is  a  symmetric  negative  definite  matrix.  (See 
Appendix  F  for  a  proof  that  Q'  is  negative  definite.)  Note  that  E(x)  =  E'  (x)  for  all 
{Q 1}  • 

The  global  minimum  of  a  strictly  concave  quadratic  function  is  attained  at  an 
extreme  point  of  the  feasible  polytope.  By  relaxing  the  integer  restrictions  in  the  polytope 
P  above  to  get  the  new  polytope  P'  consisting  of  the  constraints  (Equation  8.1)  and  (Equa¬ 
tion  8.2),  and  the  bounds  0<x{j<  1  for  i=l,...,n  and  it  can  be  seen  that  the 

extreme  points  of  P'  correspond  to  the  feasible  points  of  P.  Hence,  a  global  minimum  for 
the  strictly  concave  quadratic  problem: 


min  E'  (x) 

X  €  P' 


Equation  8.8 


will  also  be  a  global  minimum  of  the  discrete  quadratic  assignment  problem  (Equation 

8.6). 

Thus,  the  original  quadratic  assignment  problem  is  equivalent  to  the  minimization 
of  a  concave  quadratic  function  on  the  unit  hypercube  with  nN  variables,  each  of  which  is 
restricted  to  the  interval  [0,1].  Two  different  computational  methods  for  this  class  of  prob¬ 
lem  have  recently  oeen  developed.  Either  of  these  methods,  linear  underestimating  func¬ 
tions  or  a  stochastic  method  solving  multiple  cost  row  linear  programs,  can  be  used  to 
provide  the  minimum  of  the  concave,  quadratic  function.  The  global  or  near  global  solu¬ 
tion  to  »his  problem  provides  a  convenient  starting  point  for  the  relaxed,  local  minimiza¬ 
tion  problem. 
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The  solution  of  the  quadratic  assignment  problem  provides  the  conformation  with 
the  lowest  potential  energy  over  all  possible  conformations  on  the  lattice.  Since  this  con¬ 
formation  may  not  lie  at  the  exact  position  of  the  true  global  minimum,  a  local  minimiza¬ 
tion  procedure  can  be  used  to  find  the  true  solution.  This  minimization  is  therefore  an 
unconstrained  local  minimization  with  3(n-l)  variables  (Recall  that  are 

fixed.)  over  a  continuous  potential  energy  function  (See  Chapter  II.).  Any  of  the  local  min¬ 
imization  methods  described  in  Chapter  VI,  including  the  quasi-Newton  method,  will  sat¬ 
isfy  the  requirements  of  this  problem. 


D.  Potential  for  Lattice  Based  Solutions 

Unfortunately,  any  minimization  based  on  the  quadratic  assignment  method  is 

totally  reliant  on  the  selection  of  an  appropriate  lattice.  The  establishment  of  a  lattice 
which  can  accurately  model  the  system  quickly  forces  the  problem  to  become  intractable 
as  problem  size  increases.  Even  for  a  very  small  problem,  the  r  astruction  of  a  lattice 
which  can  accurately  represent  the  molecule  in  three-dimensional  space  incurs  a  very 
large  number  of  variables  (See  Table  8.1.).  Computing  a  solution  based  on  such  a  large 

Table  8.1  Number  of  Variables  Required  for  Quadratic  Assignment  Formulation 


Lattice  Size 


Number  of 
Beads 

5x5x5 

6x6x6 

7x7x7 

10x10x10 

5 

625 

1080 

1715 

5000 

6 

750 

1296 

2058 

6000 

7 

875 

1512 

2401 

7000 

12 


1500 


2592 


4116 


12000 
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number  of  variables  is  currently  intractable.  Even  the  simplest  problem  described  in  Table 
8.1  (5  beads  on  a  5x5x5  lattice),  required  131  minutes  to  produce  a  minimum  energy  con¬ 
figuration  on  the  lattice.  Once  this  solution  is  given,  the  local  minimization  phase  must 
still  be  run,  and  the  final  configuration  remains  heavily  dependent  on  the  selection  of  the 
lattice.  As  problem  size  increases,  the  number  of  variables  required  to  form  a  suitable  lat¬ 
tice  increases  at  a  rate  which  can  be  described  by  0(n4).  The  625  variable  problem 
incurred  by  the  5  bead  hydrocarbon  already  approaches  the  limits  of  current  optimization 
methods.  Without  prior  knowledge  of  the  protein  structure  (upon  which  an  appropriate  lat¬ 
tice  could  be  based),  the  quadratic  assignment  formulation  will  not  (Phillips  and  Rosen 
1994)  provide  satisfactory  results. 
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IX.  Simulated  Annealing 

A.  Optimization  Analogous  to  Physical  Annealing 

Optimization  by  simulated  annealing  is  one  of  the  newest  and  most  promising 

methods  of  finding  the  minimum  or  maximum  of  a  function  with  multiple  local  optima. 
Simulated  annealing  is  analogous  to  the  physical  process  of  annealing  in  condensed  matter 
physics.  In  physical  annealing,  a  solid  is  heated  until  it  melts.  It  is  then  carefully  cooled 
until  it  crystallizes  into  a  state  with  a  perfect  lattice.  During  this  process  the  free  energy  of 
the  solid  is  minimized.  Likewise,  simulated  annealing  attempts  to  optimize  a  function  by 
reducing  the  “temperature”  of  the  solution’s  potential  function  until  a  global  minimum 
energy  is  reached  (Laarhoven  and  Aarts  1989). 

B.  Simulated  Annealing  Theory 

The  name  “simulated  annealing”  originated  from  the  analogy  with  the  physical 
annealing  process  of  solids.  Annealing  consists  of  two  steps  (Aarts  and  Korst  1989): 

•  Increase  the  temperature  to  a  maximum  value  at  which  time  the  solid 

melts. 

•  Decrease  carefully  the  temperature  until  the  particles  arrange  themselves 

into  the  ground  state  of  the  solid. 

The  solid  shifts  from  a  state  in  which  all  particles  arrange  themselves  randomly,  to  a 
ground  state  in  which  the  particles  are  arranged  in  a  highly  structured  lattice.  In  this  lattice 
state,  the  energy  of  the  system  is  minimal.  Thus,  the  ground  state  is  obtained  only  if  the 
maximum  temperature  is  sufficiently  high,  and  the  cooling  is  done  sufficiently  slow. 

A  computational  equivalent  of  simulated  annealing  can  be  based  on  the  following 


ideas: 


Solutions  in  a  combinatorial  optimization  problem  are  equivalent  to 
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states  in  a  physical  system. 

•  The  cost  of  a  solution  is  equivalent  to  the  energy  of  a  state. 

The  role  of  temperature  is  represented  by  a  control  parameter,  c.  If  i  and  j  represent  two 

solutions  to  a  combinatorial  optimization  problem,  with  corresponding  costs  f(i),  and  f (j), 

then  the  probability  that  solution  j  is  accepted  after  visiting  solution  i  is: 

1  if  / O')  */(/) 

(. m-m ) 

*  c  if /O')  >/(/)  • 

For  every  solution  generated,  this  acceptance  criterion  is  applied,  giving  simulated  anneal¬ 
ing  the  unique  feature  of  being  able  to  accept  limited  deteriorations  in  potential  value.  The 
simulated  annealing  algorithm  can  be  stated  as: 

1 .  Given  some  initial  configuration,  <J>q,  with  corresponding  control  parameter, 

c0,  and  number  of  transitions,  Lk  and  set  it  4-  0  and  / (0)  =  F  (0O)  . 

2.  Set  t  <—  1 . 

3.  Generate  a  new  point  <|>*,  with  potential,  f(<j)*). 

4.  If  /(<}>*)  </($)  then  <t>  = 

5.  Otherwise,  if  e  is  greater  than  a  random  number  on  the  interval 

[0,1)  then  <>  =  <}>*. 

6.  Set  1 4-  /  +  1  and  if  t  is  less  than  L*  go  back  to  step  2. 

7.  Calculate  new  number  of  transitions,  Lk. 

8.  Calculate  new  control  parameter,  ck. 

9.  Return  to  step  2. 

In  this  algorithm,  ck  is  the  value  of  the  control  parameter,  and  L*  the  number  of  transitions 
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generated  at  the  iteration  of  the  algorithm  (Aarts  and  Korst  1989).  The  generation  of 
the  next  set  of  parameters,  represented  by  step  three  above  is  done  by  applying  a  perturba¬ 
tion  mechanism.  This  mechanism  transforms  the  current  state  into  the  next  state  by  a  small 
distortion.  In  the  analogy  with  the  physical  process  of  annealing,  this  change  may  be  the 
displacement  of  a  particle.  In  combinatorial  optimization,  this  change  represents  slight 
parameter  alterations. 

Simulated  annealing  is  essentially  a  generalization  of  a  local  search.  When  the 
control  parameter  is  set  to  zero,  simulated  annealing  is  identical  to  a  local  search  algo¬ 
rithm,  since  the  value  of  the  potential  function  is  only  allowed  to  decrease.  However,  the 
ability  to  move  through  ranges  of  higher  potential  function  values  enables  simulated 
annealing  to  find  high  quality  solutions  which  do  not  strongly  depend  on  the  initial  solu¬ 
tion.  Thus,  simulated  annealing  has  the  capability  to: 

•  process  potential  functions  possessing  arbitrary  degrees  of  nonlinearities, 

discontinuities,  and  stochasticity. 

•  process  arbitrary  boundary  conditions  and  constraints  imposed  on  these 

potential  functions. 

•  be  implemented  with  a  minimal  degree  of  coding  relative  to  other  nonlin¬ 

ear  optimization  algorithms. 

•  statistically  guarantee  the  finding  of  an  optimal  solution  (ergodicity). 

(Ingber  1993) 

C.  Simulated  Annealing  and  the  Protein  Folding  Problem 

Implementations  of  simulated  annealing  on  the  protein  folding  problem  did  not 

meet  expected  results.  An  implementation  of  the  simulated  annealing  algorithm  on  the 
trivially  simple  4  bead  hydrocarbon,  butane  (containing  three  local  minima)  failed  to 
return  the  global  minimum  after  several  thousand  function  evaluations.  Simulated  anneal¬ 
ing  requires  the  accurate  selection  of  temperature  schedules,  control  parameters,  and  stop- 
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ping  criteria.  With  no  knowledge  of  the  potential  function  structure  or  complexity, 
selection  of  these  values  becomes  extremely  difficult.  To  further  complicate  the  situation, 
the  user  must  provide  the  perturbation  mechanism  used  to  generate  new  points.  Even 
extremely  simple  problems,  such  as  a  concave  quadratic,  become  very  slow  if  simulated 
annealing  is  not  set  up  correctly.  For  instance,  even  if  the  global  minimum  was  selected  on 
the  first  attempt,  the  simulated  annealing  algorithm  will  still  spend  significant  amounts  of 
time  “bouncing”  around  the  minimum  point,  occasionally  allowing  increases  in  function 
value,  only  to  insure  that  the  true  global  minimum  has  been  found.  With  correct  parameter 
choices,  the  simulated  annealing  algorithm  may  be  able  to  make  this  decision  quickly,  but 
without  prior  knowledge  of  the  function,  simulated  annealing’s  sensitivity  to  parameter 
changes  makes  its  use  prohibitive. 


X.  Conclusions 


Three  computational  methods  for  solving  the  protein  folding  problem  were  exam¬ 
ined.  Two  of  these  methods  failed  to  produce  acceptable  results.  A  lattice  based  method, 
which  allows  the  formulation  of  a  quadratic  assignment  problem,  failed  due  to  the 
immense  number  of  variables  required  to  accurately  model  a  protein  chain  of  more  than 
trivial  length  (more  than  five  amino  acids).  Solving  such  a  problem,  containing  a  very 
large  number  of  variables,  is  currently  computationally  intractable. 

Simulated  annealing  has  the  theoretical  ability  to  generate  solutions  to  the  prob¬ 
lem,  but  is  completely  dependent  on  the  selection  of  parameters  which  affect  its  operation. 
Without  knowledge  of  the  characteristics  of  the  potential  energy  function,  the  selection  of 
these  values  is  very  difficult. 

The  only  successful  method  tested  was  the  convex  quadratic  global  underestima¬ 
tor.  Using  this  technique,  predictions  for  n-chain  hydrocarbons  of  length  5  through  15  and 
22  were  made.  However,  large  increases  in  local  minimization  times  make  predictions  for 

chains  of  length  greater  than  22  unlikely  using  the  current  technique.  The  molecules  tested 

3  2 

produced  a  least  squares  fit  of  0.23n  -  4.88n  +  36.03n  -  85.8836  for  total  prediction  time 
(minutes).  Using  this  approximating  function,  computation  of  an  protein  sized  molecule 
(100  residues),  would  take  133  days  of  CPU  time.  For  the  underestimator  to  realize  its 
potential,  a  superior  method  of  finding  local  minima  must  be  developed.  The  local  mini¬ 
mization  procedures  rely  on  many  function  evaluations.  Improvements  in  the  function 
evaluation  time  (including  conversions)  would  be  of  immense  benefit,  and  perhaps  pro¬ 
vide  the  key  to  reducing  local  minimization  times.  Recent  work  by  Greengard  (1988)  has 
attempted  to  specifically  address  this  problem. 
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Initial  Position,  y 


Minimized  Position,  <)) 


180.116728 


-0.215524 


179.627701  180.119720  li 


179.881415  180.000000 


-0.215452 


179.630935  -61689425 


180.000000 


-0.215452 


180.407944  61646775 


179.718991  -61709134  178.586693  -61896149 


-0.195689 


178.613080  -61715267  179.718298 


62.919639 


62.70792S  180.299240 


61908155  -178.598808  61670638 


180.321668 


61895050 


0.195526 


62.500992  180.750770 


-0.139428 


-62.475940  179159740  181740408  62.475995 


-0.139423 


64.173448  63.997801  180.864636  179.735869 


XII.  Appendices 


A.  Conformers  for  Heptane 


Table  12.1  Local  Minima  from  Preferred  Angles,  Heptane 


•0.123230 


180 

180 

180 

180 

180 

-60 

180 

180 

180 

180 

-60 

180 

180 

180 

60 

180 

180 

60 

180 
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180 

60 

60 
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60 
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60 
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60 

180 

60 

-60 

60 

-60 

180 

60 

180 

60 

60 

180 

180 

-60 

-60 

180 

180 

60 

60 

60 

180 

180 

•179.715697 
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Table  12.1  Local  Minima  from  Preferred  Angles,  Heptane 


U«t» 

Initial  Position,  <j>^ 

Minimized  Position,  <J>* 

— 

-0.125181 

-60 

60 

-60 

-60 

-179.762381 

179.114165 

-63.985617 

-64.185171 

-0.125144 

180 

180 

60 

60 

179.760412 

180.856949 

63.972239 

64.164748 

-0.125125 

-60 

-60 

180 

180 

-64.177304 

-64.007009 

179.159146 

180.248780 

-0.125125 

-60 

-60 

60 

-60 

-64.210703 

-64.018941 

179.154209 

-179.752857 

-0.125101 

180 

180 

-60 

-60 

180.222763 

179.150008 

-64.035518 

-64.156472 

-0.124850 

-60 

180 

60 

180 

-62.737240 

180.041359 

62.538518 

180.313344 

-0.124826 

60 

180 

-60 

180 

62.766715 

179.950797 

-62.544082 

179.694691 

-0.124723 

180 

-60 

180 

60 

179.676542 

-62.513659 

179.959773 

62.721797 

-0.124681 

180 

60 

180 

-60 

180.321851 

62.542211 

180.033777 

-62.749200 

-0.115072 

-60 

-60 

180 

-60 

-63.801675 

-64.208098 

177.402217 

-62.874698 

-0.115043 

60 

180 

60 

60 

62.902069 

182.598851 

64.182034 

63.817398 

-0.115031 

60 

60 

180 

60 

63.819117 

64.207828 

182.636566 

62.861170 

-0.114919 

-60 

180 

-60 

-60 

-62.897443 

177.358042 

-64.224387 

-63.846126 

-0.104232 

60 

180 

180 

60 

62.576398 

180.185008 

180.185012 

62.575741 

-0.104125 

-60 

180 

180 

-60 

-62.534582 

179.808088 

179.807766 

-62.603440 

-0.104118 

60 

60 

-60 

60 

62.550498 

180.191468 

-179.798849 

62.592457 

0.009673 

180 

60 

60 

60 

180.725871 

63.726114 

64.947038 

64.047553 

0.009674 

60 

60 

60 

180 

64.016389 

64.958687 

63.664674 

180.702491 

0.009691 

180 

-60 

-60 

-60 

179.273906 

-63.686327 

-64.964465 

-64.041972 

0.009736 

-60 

-60 

-60 

180 

-64.055443 

-64.987973 

-63.684967 

179.269252 

0.045243 

60 

180 

-60 

-60 

62.685626 

179.523486 

-64.182031 

-64.199601 

0.045296 

-60 

-60 

180 

60 

-64.176786 

-64.166840 

179.514057 

62.669550 

0.045313 

-60 

180 

60 

60 

-62.666795 

180.526585 

64.165301 

64.223136 

0.045388 

60 

60 

180 

-60 

64.209791 

64.193968 

180.526278 

-62.685254 

0,249514 

60 

60 

60 

60 

64.088236 

64.918552 

64.907186 

64.083887 

0.249592 

-60 

-60 

-60 

-60 

-64.088099 

-64.893060 

-64.880987 

-64.070387 

63 


Table  12.1  Local  Minima  from  Preferred  Angles,  Heptane 


U«t» 

Initial  Position,  <J>^ 

Minimized  Position,  <{>* 

2.026672 

60 

60 

-60 

-60 

63.000643 

79.493002 

-79.562918 

-62.995198 

2.026839 

-60 

-60 

60 

60 

-62.996259 

-79.546417 

79.567594 

62.981162 

2.058274 

-60 

60 

60 

180 

•80.186256 

78.462644 

62.066901 

181.075020 

2.058303 

60 

-60 

-60 

180 

80.184791 

-78.395394 

-62.095009 

178.921228 

1058343 

180 

60 

60 

-60 

181.106215 

62.088491 

78.404974 

-80 

.180360 

2.058378 

180 

-60 

-60 

60 

178.917917 

-62.0S23Q5 

-78.445940 

80 

.186494 

2.070151 

180 

-60 

60 

180 

179.287701 

-79.399090 

79.614853 

180.690834 

1070219 

180 

60 

-60 

180 

180.714220 

79.494318 

-79365457 

179.293132 

1097911 

180 

60 

-60 

-60 

180.595288 

79.299839 

-79.307316 

-62.430627 

2.097949 

-60 

-60 

60 

180 

-62.396732 

-79.336814 

79.286330 

180.623456 

2.097992 

180 

-60 

60 

60 

179.378780 

-79.278332 

79.236470 

62.423103 

1098034 

60 

60 

-60 

180 

62411331 

79.279927 

-79.366127 

179.373374 

2.158417 

180 

180 

-60 

60 

180.225655 

179.209800 

-79.355190 

79.491831 

2.158423 

-60 

60 

180 

180 

-79.398437 

79.611536 

180.835907 

179.812488 

2.158439 

180 

180 

60 

-60 

179816175 

180.791002 

79.434474 

-79.543133 

2.158467 

60 

-60 

180 

180 

79.514152 

-79.455677 

179.243947 

180.205755 

2.230068 

60 

180 

60 

-60 

62.935288 

182.279514 

82.013853 

-77.339439 

2.230082 

-60 

60 

180 

60 

-77.265866 

82.102048 

182.245218 

62.944156 

2.230124 

60 

-60 

180 

-60 

77.368066 

-81.914342 

177.749227 

-62.976229 

2.230226 

-60 

180 

-60 

60 

-62.947196 

177.778962 

-81.899852 

77.352625 

2.336922 

60 

60 

60 

-60 

64.085833 

63.375309 

78.562093 

-79.969246 

2.337016 

-60 

60 

60 

60 

-79.963708 

. 

78370336 

63.378969 

64.079620 

2.337018 

60 

-60 

-60 

-60 

79.875053 

-78.515297 

-63.341112 

-64.066549 

2.337127 

-60 

-60 

-60 

60 

-64.098840 

-63.359511 

-78.518880 

79.935753 

2.343269 

-60 

60 

180 

-60 

-79.299852 

79.719299 

180.326332 

-62.698935 

2.343326 

60 

180 

-60 

60 

62.700984 

179.692084 

-79.670324 

79.291646 
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Table  12.1  Local  Minima  from  Preferred  Angles,  Heptane 


u«t>) 

Initial  Position,  <f>^ 

Minimized  Position,  <j>* 

2.343330 

60 

-60 

180 

60 

79.431571 

-79.717732 

179.674046 

62.754099 

2.343345 

-60 

180 

60 

-60 

-62.719402 

180.352229 

79.647337 

-79.369962 

4.151268 

60 

-60 

-60 

60 

78.132822 

-76.957224 

-76.943581 

78.178186 

4.151353 

-60 

60 

60 

-60 

-78.142142 

76.959901 

76.950608 

-78.166430 

4.384072 

180 

-60 

60 

-60 

178.741629 

-93.723491 

69.284953 

-94.158124 

4.384075 

-60 

60 

-60 

180 

-94.198578 

69.340226 

-93.775408 

178.766022 

4.384250 

60 

-60 

60 

180 

94.299939 

-69.383118 

93.736766 

181.247700 

4.502235 

60 

-60 

60 

60 

95.723592 

-70.556034 

90.568777 

61.256347 

B.  Global  Underestimator  Run:  Butane  (4  beads) 

Host  Machine:  haney  Tue  Apr  12  10:46:53  EDT  1994 

How  many  randomizations  (at  least  3):  10 

Within  what  radius  of  phi[l]  =  180.0000  (degs):  180 

— >  Using  10  randomizations 

running  QN-BFGS  (using  all  local  minima) ...  pass  #1 

Elapsed  Real  Time  for  QN-BFGS  =  1 .000  Seconds 

— >  1 1  points  reduced  to  1 1  unique  conformers 

solving  LP ... 


Conformations:  11  (possibly  non-unique) 

Dihedral  Angles:  I 

Quad  Fit  Error  =  0,014714  kcal/mol  (sum  of  errors) 
LP  Pivots  =  4 

Interpolated  Conformations: 


I  2  3  (3  out  of  1 1 ) 
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Convex  Quadratic  Underestimator: 
c[0]  =  +0.553564 

c[l]  = -0.357722  d[l]  =  +0.113864 

Max  error  in  Quad  Fit  =  0.006454  kcal/mol 

Distribution  of  Errors  (%  of  Max) 

10%  20%  30%  40%  50%  60%  70%  80%  90%  100% 

8  0  0  1  0  0  0  0  1  1 

Ranges  on  the  Dihedrals  (Degrees): 

63.01 1900  <=  phi[  1  ]  <=  297.081800 

Elapsed  Real  Time  for  LP  =  0.000  Seconds 

LP  Predicted  Dihedrals  Current  Global  Minimum 

phi[l]  =  3.141660  rads  (180.003872  degs)  phi[l]  =  3.141593  rads  (180.000000  degs) 

Energy  =  -0.008356  kcal/mol 

LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  0.000068  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  0.000000  kcal/mol 

Repeat  (y/n)?  y 

How  many  randomizations:  5 

What  fraction  of  the  distance  to  the  global  min  should  be  used  [0,1]:  0 
Use  how  many  points  from  previous  run  (global  only/all):  g 

Within  what  radius  of  phi[l]  -  180.0039  (degs):  5 

~>  Using  5  randomizations 

running  QN-BFGS  (using  only  global  minimum) ...  pass  #2 

Elapsed  Real  Time  for  QN-BFGS  =  1 .000  Seconds 
~>  7  points  reduced  to  5  unique  conformers 

solving  LP  ... 


Conformations:  5  (possibly  non-unique) 
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Dihedral  Angles:  1 

Quad  Fit  Error  =  0.000681  kcal/mol  (sum  of  errors) 

LP  Pivots  =  6 

Interpolated  Conformations: 

3  4  (2  out  of  5) 

Convex  Quadratic  Underestimator: 
c[0]  = +0.221661 

c[  1  ]  =  -0.073250  d[  1  ]  =  +0.000000 

Max  error  in  Quad  Fit  =  0.000348  kcal/mol 

Distribution  of  Errors  (%  of  Max) 

10%  20%  30%  40%  50%  60%  70%  80%  90%  100% 

2  0  0  1  0  0  10  0  1 

Ranges  on  the  Dihedrals  (Degrees): 

179.196100  <=phi[l]<=  180.156200 

Elapsed  Real  Time  for  LP  =  0.000  Seconds 

LP  Predicted  Dihedrals  Current  Global  Minimum 

phi[l]  =  3.141593  rads*  (1 80.000000 degs)  phi[l]  =  3.141593  rads  ( 1 80.000000  degs) 

Energy  =  -0.008356  kcal/mol 

LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  0.000000  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  0.000000  kcal/mol 


Successive  LP  predicted  dihedrals  change  (inf-norm)  =  0.003900  degs 
Successive  global  conformers  change  (inf-norm)  =  0.000000  degs 
Eneigy  improvement  =  0.000000  kcal/mol 

Repeat  (y/n)?  n 
Run  Complete.... 

C.  Global  Underestimator  Run:  Heptane  (7  beads) 


Host  Machine:  haney 


Tue  Apr  12  11:18:57  EDT  1994 
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How  many  randomizations  (at  least  9):  25 

Within  what  radius  of  phi[  1  ]  =  20.0000  (degs):  1 80 
Within  what  radius  of  phi[2]  =  330.0000  (degs):  1 80 
Within  what  radius  of  phi[3]  =  278.0000  (degs):  1 80 
Within  what  radius  of  phi[4]  =  -60.0000  (degs):  180 

-->  Using  25  randomizations 

running  QN-BFGS  (using  all  local  minima) ...  pass  #1 


Elapsed  Real  Time  for  QN-BFGS  =  40.000  Seconds 
-->  26  points  reduced  to  26  unique  conformers 


solving  LP ... 


Conformations:  26  (possibly  non-unique) 

Dihedral  Angles:  4 

Quad  Fit  Error  =  24.465463  kcal/mol  (sum  of  errors) 

LP  Pivots  =  13 

Interpolated  Conformations: 

1  8  9  10  12  13  17  21  26  (9  out  of  26) 

Convex  Quadratic  Underestimate* : 

c[0]  = +1.313334 

c[  1  ]  =  -0.299 179  d[l]  =  +0.087703 

c[2]  =  -0.273989  d[2]  =  +0.083520 

c[3]  =  -0.225441  d[3]  =  +0.079214 

c[4]  =  -0.344606  d[4]  =  +0. 100744 

Max  error  in  Quad  Fit  =  4.677598  kcal/mol 

Distribution  of  Errors  (%  of  Max) 

10%  20%  30%  40%  50%  60%  70%  80%  90% 

17  00070000 

Ranges  on  the  Dihedrals  (Degrees): 

62.664800  <=  phi[l]  <=  297.534500 
62.362800  <=  phi{2]  <=  297.492700 
62.576500  <=  phi[3]  <=  297.954800 
62.516700  <=  phi[4]  <=  297.321200 


100% 

2 


Elapsed  Real  Time  for  LP  =  1 .000  Seconds 


LP  Predicted  Dihedrals  Current  Global  Minimum 


phi[  1  ]  =  3.411274  rads  (195.451627  degs) 
phi[2]  =  3.280520  rads  (187.959930  degs) 
phi[3]  =  2.845974  rads  (163.062310  degs) 
phi[4]  =  3.42061 1  rads  (195.986554  degs) 


phi[l]  =  3.143848  rads  (180.129200  degs) 
phi[2]  =  3.133486  rads  ( 1 79.535500  degs) 
phi[3]  =  5.192398  rads  (297.502500  degs) 
phi[4]  =  3.135182  rads  (179.632700  degs) 

Energy  =  -0.296754  kcal/mol 


LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  2.346424  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  0.226207  kcal/mol 


Repeat  (y/n)?  y 

How  many  randomizations:  20 

What  fraction  of  the  distance  to  the  global  min  should  be  used  [0,1  ]:  0 
Use  how  many  points  from  previous  run  (global  only/all):  g 

Within  what  radius  of  phi [  1  ]  =  195.4516  (degs):  20 
Within  what  radius  of  phi[2]  =  187.9599  (degs):  15 
Within  what  radius  of  phi[3]  =  163.0623  (degs):  140 
Within  what  radius  of  phi[4]  =  195.9866  (degs):  20 

-->  Using  20  randomizations 

running  QN-BFGS  (using  only  global  minimum) ...  pass  #2 


Elapsed  Real  Time  for  QN-BFGS  =  15.000  Seconds 
— >  22  points  reduced  to  22  unique  conformers 

solving  LP ... 


Conformations:  22  (possibly  non-unique) 

Dihedral  Angles:  4 

Quad  Fit  Error  =  0.01 1400  kcal/mol  (sum  of  errors) 
LP  Pivots  =  20 

Interpolated  Conformations: 

1  2  78  11  20  21  (7  out  of  22) 

Convex  Quadratic  Underestimator: 


c[0]  =  +46.835326 
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c[  1  ]  =  -0.056420  d[  1  ]  =  +0.000000 

c[2]  =  -0.05 1 352  d[2]  =  +0.000000 

c[3]  =  -0.063887  d[3]  =  +0.0203 1 9 

c[4]  =  -29.772245  d[4]  =  +9.482895 

Max  error  in  Quad  Fit  =  0.002471  kcal/moi 

Distribution  of  Errors  (%  of  Max) 

10%  20%  30%  40%  50%  60%  70%  80%  90%  100% 

10  3223  10001 

Ranges  on  the  Dihedrals  (Degrees): 

179.383400  <=  phi[l]  <=  180.446000 
179.018300  <=  phi[2]  <=  180.749100 
62.158100  <=  phi[3]  <=  297.544400 
179.528500  <=  phi[4]  <=  181.045700 


Elapsed  Real  Time  for  LP  =  1 .000  Seconds 


LP  Predicted  Dihedrals  Current  Global  Minimum 

phi[  1  ]  =  3.141902  rads*  (180.017700  degs)  phi[l]  =  3.141902  rads  (180.017700  degs) 

phi[2]  =  3.141477  rads*  (179.993400  degs)  phi[2]  =  3.141477  rads  (179.993400  degs) 

phi[3]  =  3.144200  rads  (180.149391  degs)  phi[3]  =  3.142771  rads  (180.067500  degs) 

phi[4]  =  3.139573  rads  (179.884306  degs)  phi[4]  =  3.141183  rads  (179.976500  degs) 

Energy  =  -0.339768  kcal/mol 

LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  0.001609  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  0.000012  kcal/mol 


Successive  LP  predicted  dihedrals  change  (inf-norm)  =  17.087100  degs 
Successive  global  conformers  change  (inf-norm)  =  1 17.435000  degs 
Energy  improvement  =  0.043014  kcal/mol 


Repeat  (y/n)?  y 

How  many  randomizations:  10 

What  fraction  of  the  distance  to  the  global  min  should  be  used  [0,1]:  0 
Use  how  many  points  from  previous  run  (global  only/all):  g 

Within  what  radius  of  phifl]  =  180.0177  (degs):  4 
Within  what  radius  of  phi[2]  =  179.9934  (degs):  4 
Within  what  radius  of  phi[3]  =  180.1494  (degs):  5 


Within  what  radius  of  phi[4]  =  179.8843  (degs):  5 
-->  Using  10  randomizations 

running  QN-BFGS  (using  only  global  minimum) ...  pass  #3 
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Elapsed  Real  Time  for  QN-BFGS  =  3.000  Seconds 
— >  12  points  reduced  to  12  unique  conformers 


solving  LP ... 


Conformations:  12  (possibly  non-unique) 

Dihedral  Angles:4 

Quad  Fit  Error  =0.000950  kcal/mol  (sum  of  errors) 

LP  Pivots  =  10 

Interpolated  Conformations: 

2  3  47  10  12  (6  out  of  12) 

Convex  Quadratic  Underestimator: 

c[0]  =  +88.986822 

c[  1  ]  =  -0.046979  d[  1  ]  =  +0.000000 

c[2]  =  -56.991488  d[2]  =  +18.152711 

c[3]  =  +0.006854  d[3]  =  +0.000000 

c[4]  =  +0.083894  d[4]  =  +0.000000 

Max  error  in  Quad  Fit  =  0.000436  kcal/mol 

Distribution  of  Errors  (%  of  Max) 

10%  20%  30%  40%  50%  60%  70%  80%  90% 

7  0  1  2  0  10  0  0 

Ranges  on  the  Dihedrals  (Degrees): 

179.194700  <=phi[l]<=  180.700300 
178.951300  <=  phi(2]  <=  180.893700 
179.722400  <=  phi[3]  <=  180.467900 
179.851200  <=  phi[4]  <=  180.888000 


Elapsed  Real  Time  for  LP  =  0.000  Seconds 


100% 

1 


LP  Predicted  Dihedrals 


Current  Global  Minimum 
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phi[l]  =  3.141902  rads*  (180.017700  degs) 
phi[2]  =  3.139558  rads  (179.883420  degs) 
phi[3]  =  3.143699  rads*  (180.120700  degs) 
phi[4]  =  3.139963  rads*  (179.906600  degs) 


phi[  1  ]  =  3.141902  rads  (180.017700  degs) 
phi[2]  =  3.141477  rads  (179.993400  degs) 
phi[3]  =  3. 143699  rads  ( 1 80. 1 20700  degs) 
phi [4]  =  3.139963  rads  (179.906600  degs) 

Energy  =  -0.339793  kcal/mol 


LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  0.001920  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  0.000033  kcal/mol 


Successive  LP  predicted  dihedrals  change  (inf-norm)  =  0.1 10000  degs 
Successive  global  conformers  change  (inf-norm)  =  0.069900  degs 
Energy  improvement  =  0.000025  kcal/mol 


Repeat  (y/n)?  y 

How  many  randomizations:  5 

What  fraction  of  the  distance  to  the  global  min  should  be  used  [0,1]:  0 
Use  how  many  points  from  previous  run  (global  only/all):  g 

Within  what  radius  of  phi[l]  =  180.0177  (degs):  0 
Within  what  radius  of  phi[2]  =  1 79.8834  (degs):  5 
Within  what  radius  of  phi[3J  =  180.1207  (degs):  0 
Within  what  radius  of  phi[4]  =  179.9066  (degs):  0 

-->  Using  5  randomizations 

running  QN-BFGS  (using  only  global  minimum) ...  pass  #4 


Elapsed  Real  Time  for  QN-BFGS  =  1 .000  Seconds 
— >  7  points  reduced  to  7  unique  conformers 


solving  LP ... 


Conformations:  7  (possibly  non-unique) 

Dihedral  Angles:  4 

Quad  Fit  Error  =  0.001902  kcal/mol  (sum  of  errors) 
LP  Pivots  =  12 

Interpolated  Conformations: 

12  347  (5  out  of  7) 


Convex  Quadratic  Underestimator: 
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c[OJ  =  +6.716253 

c[l]  =  -0.830429  d[l]  =  +0.000000 

c(2]  =  -0.06 1182  d[2]  =  +0.000000 

c[3]  =  -0.464009  d[3]  =  +0.000000 

c[4]  =  -0.890460  d[4]  =  +0.000000 

Max  error  in  Quad  Fit  =  0.001696  kcal/mol 

Distribution  of  Errors  (%  of  Max) 

10%  20%  30%  40%  50%  60%  70%  80%  90%  100% 

5  100000001 

Ranges  on  the  Dihedrals  (Degrees): 

1 79.979400  <=  phi[  1  ]  <=  180.017700 
179.261800  <=  phi[2]  <=  180.713700 
179.971500  <=phi[3]<=  180.120700 
179.906600  <=  phi [4]  <=  180.009200 


Elapsed  Real  Time  for  LP  =  0.000  Seconds 


LP  Predicted  Dihedrals  Current  Global  Minimum 

phi [  1  ]  =  3.141902  rads*  (180.017700 degs)  phi[l]  =  3.141902  rads  (180.017700  degs) 

phi[2]  =  3.141477  rads*  (179.993400  degs)  phi[2]  =  3.141477  rads  (179.993400  degs) 

phi[3J  =  3.143699  rads*  (180. 120700  degs)  phi[3]  =  3.143699  rads  (180. 120700  degs) 

phi[4]  =  3.139963  rads*  (179.906600  degs)  phi[4]  =  3. 139963  rads  (179.906600  degs) 

Energy  =  -0.339793  kcal/mol 

LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  0.000000  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  0.000000  kcal/mol 


Successive  LP  predicted  dihedrals  change  (inf-norm)  =  0.1 10000  degs 
Successive  global  conformers  change  (inf-norm)  =  0.000000  degs 
Energy  improvement  =  0.000000  kcal/mol 


Repeat  (y/n)?  n 
Run  Complete.... 

D.  Global  Underestimator  Run:  12  beads 


Host  Machine,  csserverc 


Fri  Apr  15  08:43:35  EDT  1994 
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How  many  randomizations  (at  least  19):  60 


Within  what  radius  of  phifl]  = 
Within  what  radius  of  phi[2]  = 
Within  what  radius  of  phi[3]  = 
Within  what  radius  of  phi[4]  = 
Within  what  radius  of  phi[5]  = 
Within  what  radius  of  phi[6]  = 
Within  what  radius  of  phi[7]  = 
Within  what  radius  of  phi[8]  = 
Within  what  radius  of  phi[9]  = 

— >  Using  60  randomizations 


330.0000  (degs):  180 
278.0000  (degs):  180 
123.0000  (degs):  180 
31.0000  (degs):  180 
160.0000  (degs):  180 
97.0000  (degs):  180 
27.0000  (degs):  180 
232.0000  (degs):  180 
38.0000  (degs):  180 


running  QN-BFGS  (using  all  local  minima) ...  pass  #1 


Elapsed  Real  Time  for  QN-BFGS  =  801.000  Seconds 
-->  61  points  reduced  to  61  unique  conformers 


solving  LP ... 


Conformations:  61  (possibly  non-unique) 

Dihedral  Angles:  9 

Quad  Fit  Error  =  84.557709  kcal/mol  (sum  of  errors) 
LP  Pivots  =  50 


Interpolated  Conformations: 

2  5  13  20  21  23  25  29  31  33  37  39  42  55  57  61  (16  out  of  61) 


Convex  Quadratic  Underestimator: 


c[0]  =  +4.689026 
c[l]  = -0.211345 
c[2]  = -1.757095 
c(3]  = -1.230106 
c[4]  = +0.149471 
c[5]  =  -0.109763 
c[6)  =  -0.016002 
c[7]  =  +0.012777 
c[8]  =  -0.931717 
c[9]  =  +0.133952 


d[l]  =  +0.083985 
d[2]  =  +0.538594 
d[3]  =  +0.284693 
d[4]  =  +0.000000 
d[5]  =  +0.000000 
d[6)  =  +0.025305 
d[7]  =  +0.000000 
d[8)  = +0.1 88731 
d[9]  =  +0.003914 


Max  error  in  Quad  Fit  =  6.538950  kcal/mol 


Distribution  of  Errors  (%  of  Max) 


10%  20%  30%  40%  50%  60%  70%  80%  90% 


100% 
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26  14  3  7  1  6  0  1  2  1 

Ranges  on  the  Dihedrals  (Degrees): 

60.391500  <=  phi(l]  <=  297.814200 
61.145700  <=  phi[2]  <=  297.677800 
60.504700  <=  phi(3J  <=  299.337000 
60.057200  <=  phi[4]  <=  299.964500 
61.774000  <=  phi(5]  <=  301.509700 
60.360500  <=  phi(6]  <=  300.002300 
61.171200  <=  phi(7]  <=  298.399200 
59.886200  <=  phi[8]  <=  299.800300 
6 1 .052 1 00  <=  phi[9]  <=  298.280 1 00 


Elapsed  Real  Time  for  LP  =  1 .000  Seconds 


LP  Predicted  Dihedrals  Current  Global  Minimum 


phi[l]  =  2.516461  rads  (144.182610  degs) 
phi[2]  =  3.262374  rads  (186.920255  degs) 
phi[3]  =  4.320816  rads  (247.564507  degs) 
phi[4]  =  5.148753  rads*  (295.001 800  degs) 
phi[5]  =  2.881113  rads*  (165 .075600  degs) 
phi[6]  =  0.632365  rads  (  36.231854  degs) 
phi[7]  =  3.109494  rads*  (178.160900  degs) 
phi[8]  =  4.936746  rads  (282.854707  degs) 
phi[9]  =  3.475300  rads  (199.120016  degs) 


phifl]  =  5.164642  rads  (295.912200 degs) 
phi[2]  =  3.002010  rads  (172.002500  degs) 
phi[3]  =  5.187037  rads  (297.195300  degs) 
phi[4]  =  5.148753  rads  (295.001800  degs) 
phi[5]  =  2.881 1 13  rads  (165.075600  degs) 
phi[6]  =  5.153449  rads  (295.270900  degs) 
phi[7]  =  3.109494  rads  (178.160900  degs) 
phi[8]  =  3.139512  rads  (179.880800  degs) 
phi[9]  =  1.065560  rads  (  61.052100  degs) 


Energy  =  - 1 .784457  kcal/mol 


LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  4.521084  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  0.638772  kcal/mol 


Repeat  (y/n)?  y 


How  many  randomizations:  60 

What  fraction  of  the  distance  to  the  global  min  should  be  used  [0,1]:  0 
Use  how  many  points  from  previous  run  (globai  only/all):  g 


Within  what  radius  of  phi[  ‘.  ] 
Within  what  radius  of  phi[2] 
Within  what  radius  of  phi[3] 
Within  what  radius  of  phi[4] 
Within  what  radius  of  phi[5] 
Within  what  radius  of  phi[6] 
Within  what  radius  of  phi[7] 
Within  what  radius  of  phi[8] 
Within  what  radius  of  phi[9] 


144.1826  (degs):  150 
186.9203  (degs):  20 
247.5645  (degs):  60 
295.0018  (degs):  180 
165.0756  (degs):  180 
36.2319  (degs):  160 
178.1609  (degs):  180 
282.8547  (degs):  120 
199. 1200  (degs):  150 


-->  Using  60  randomizations 
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running  QN-BFGS  (using  only  global  minimum) ...  pass  #2 


Elapsed  Real  Time  for  QN-BFGS  =  765.000  Seconds 
-->  62  points  reduced  to  62  unique  conformers 

solving  LP ... 


Conformations:  62  (possibly  non-unique) 

Dihedral  Angles:  9 

Quad  Fit  Error  =  8 1 .330506  kcal/mol  (sum  of  errors) 
LP  Pivots  =  49 


Interpolated  Conformations: 

1  37  15  17  18  19  20  22  26  33  35  36  43  49  60  (16  out  of  62) 


Convex  Quadratic  Underestimator: 


c[0]  =  +1225.378595 
c[l]  =  +0.120373 
c[2]  =  -809.515465 
cl  3]  =  +0.397725 
c[4]  =  -0.388046 
c[5]  = -0.103557 
c[6]  =  -0.816778 
c[7]  = -0.821738 
c[8]  = -0.183396 
c[9]  = -1.261490 

Max  error  in  Quad  Fit 


d[l]  =  +0.000000 
d[2]  =  +267.458814 
d[3J  =  +0.000000 
d(4]  =  +0.104867 
d(5]  =  +0.000000 
d[6]  = +0.2231 16 
d[7]  =  +0.232574 
d[8]  =  +0.089597 
d[9]  =  +0.369553 

=  4.766802  kcal/mol 


Distribution  of  Errors  (%  of  Max) 


10% 

22 


20% 

11 


30% 

4 


40% 

0 


50% 

14 


60% 

5 


Ranges  on  the  Dihedrals  (Degrees): 


60.448900  <= 
172.002500  <= 
173.720700  <= 
60.647700  <= 
62.011400  <= 
59.459000  <= 
60.256700  <= 
62.973100  <= 
60.689600  <= 


phi[l]<= 

;  phi[2]  <= 
phi[3]  <= 
phi[4]  <= 
phi[5]  <= 
phi[6]  <= 
phi[7]  <= 
phi[8]  <= 
phi[9]  <= 


298.329800 

:  182.659000 

:  300.268000 

300.452200 

298.101000 

300.846400 

299.676000 

301.085700 

298.120500 


70% 

1 


80% 

0 


90% 

2 


100% 

3 
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Elapsed  Real  Time  for  LP  =  1 .000  Seconds 


LP  Predicted  Dihedrals  Current  Global  Minimum 


phil  1  ]  =  1.088288  rads*  (  62.354300  degs) 
phi[2]  =  3.026692  rads  (173.416680  degs) 
phi[3]  =  5.181839  rads*  (296.897500  degs) 
phi[4]  =  3.700363  rads  (212.015201  degs) 
phi[5]  =  2.962707  rads*  (169.750600  degs) 
phi[6]  =  3.660777  rads  (209.747092  degs) 
phi[7]  =  3.533232  rads  (202.439306  degs) 
phi[8]  =  2.046899  rads  (1 17.278668  degs) 
phi[9]  =  3.413556  rads  (195.582374  degs) 


phi[  1  ]  =  1.088288  rads  (  62.354300  degs) 
phi[2]  =  3.093175  rads  (177.225900  degs) 
phi[3]  =  5.181839  rads  (296.897500  degs) 
phi [4]  =  5.169088  rads  (296.166900  degs) 
phi[5]  =  2.962707  rads  (169.750600  degs) 
phi[6]  =  5.171863  rads  (296.325900  degs) 
phi[7]  =  5.200229  rads  (297.951200  degs) 
phi[8]  =  3.078949  rads  (176.410800  degs) 
phi[9]  =  3.161196  rads  (181.123200  degs) 


Energy  =  -2.476850  kcal/mol 


LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  1 .666997  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  1.341558  kcal/mol 


Successive  LP  predicted  dihedrals  change  (inf-norm)  =  173.515200  degs 
Successive  global  conformers  change  (inf-norm)  =  233.557900  degs 
Energy  improvement  =  0.692393  kcal/mol 


Repeat  (y/n)?  y 


How  many  randomizations:  50 

What  fraction  of  the  distance  to  the  global  min  should  be  used  [0,1]:  .5 
Use  how  many  points  from  previous  run  (global  only/all):  g 


Within 

Within 

Within 

Within 

Within 

Within 

Within 

Within 

Within 


what  radius  of  phi[  1  ] 
what  radius  of  phi[2] 
what  radius  of  phi[3] 
what  radius  of  phi[4] 
what  radius  of  phi[5] 
what  radius  of  phi[6] 
what  radius  of  phi[7] 
what  radius  of  phi[8] 
what  radius  of  phi[9] 


62.3543  (degs):  90 
175.3213  (degs):  10 
296.8975  (degs):  20 
254.0911  (degs):  50 
169.7506  (degs):  180 
253.0365  (degs):  50 
250.1953  (degs):  50 
146.8447  (degs):  50 
188.3528  (degs):  20 


>  Using  50  randomizations 


running  QN-BFGS  (using  only  global  minimum) ...  pass  #3 


Elapsed  Real  Time  for  QN-BFGS  =  526.000  Seconds 
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-->  52  points  reduced  to  52  unique  conformers 
solving  LP ... 


Conformations:  52  (possibly  non-unique) 

Dihedral  Angles:  9 

Quad  Fit  Error  =  57.460679  kcal/mol  (sum  of  errors) 
LP  Pivots  =  37 


Interpolated  Conformations: 

1  5  19  20  24  26  29  31  35  36  38  40  48  (13  out  of  52) 


Convex  Quadratic  Underestimator: 


c[0]  = +15577.752075 
c[l]  = -1.572693 
c[2]  =  -25.942274 
c[3]  =  -5963.287255 
c[4]  =  +0.145274 
c(5]  =  -5.584252 
c[6]  = +0.321134 
c[7]  =  +0.015908 
c[8]  = -1.244513 
c[9]  =  +4.039473 


d[l]  = +0.545076 
d[2]  =  +0.000000 
d[3]  = +1147.053749 
d[4]  =  +0.000000 
d{5]  = +1.658767 
d[6]  =  +0.000000 
d[7]  =  +0.000000 
d[8]  =  +0.000000 
d[9J  =  +0.000000 


Max  error  in  Quad  Fit  =  3.752148  kcal/mol 

Distribution  of  Errors  (%  of  Max) 

10%  20%  30%  40%  50%  60%  70%  80%  90% 

21  5  4  3  1  10  3  1  1 


Ranges  on  the  Dihedrals  (Degrees): 


61.855600  <= 
174.944100  <= 
295.252600  <= 
173.238200  <= 
62.177300  <= 
173.984000  <= 
173.549100  <= 
62.715800  <= 
179.081600  <= 


phi[l]<= 
phi[2]  <= 
phi[3]  <= 
phi[4]  <= 
phi[5]  <= 

:  phi [6]  <= 
:  phi[7]  <= 
phi[8]  <= 

:  Phi[9]  <= 


299.884300 
180.196500 
301.388600 
298.521300 
297.667100 
:  298.358300 
:  298.672900 
182.877300 
:  181.467600 


Elapsed  Real  Time  for  LP  =  1 .000  Seconds 


100% 

3 


LP  Predicted  Dihedrals 


Current  Global  Minimum 
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phi(  1  ]  =  2.885273  rads  (165.313959  degs) 
phi[2]  =  3.100785  rads*  (177.661900  deg) 
phi[3]  =  5.198786  rads  (297.86851 1  degs) 
phi[4]  =  5. 161232  rads*  (295.716800  degs) 
phi[5]  =  3.366508  rads  (192.886687  degs) 
phi(6]  =  5. 16051 1  rads*  (295.675500  degs) 
phi[7]  =  5.189426  rads*  (297.332200  degs) 
phi[8]  =  3.092086  rads*  (177.163500  degs) 
phi[9j  =  3.144113  rads*  (180. 144400  degs) 


phi[l]  =  3.146818  rads  (180.299400  degs) 
phi[2]  =  3.100785  rads  (177.661900  degs) 
phi[3]  =  5.191986  rads  (297.478900  degs) 
phi[4]  =  5.161232  rads  (295.716800  degs) 
phi[5]  =  2.960618  rads  (169.630900  degs) 
phi[6]  =  5.160511  rads  (295.675500  degs) 
phi[7]  =  5.189426  rads  (297.332200  degs) 
phi[8]  =  3.092086  rads  (177.163500  degs) 
phi[9]  =  3. 144 1 1 3  rads  ( 1 80. 1 44400  degs) 


Energy  =  -3.006272  kcal/mol 


LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  0.405890  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  0.181801  kcal/mol 


Successive  LP  predicted  dihedrals  change  (inf-norm)  =  102.959700  degs 


Successive  global  conformers  change  (inf-norm)  =  1 17.945100  degs 
Energy  improvement  =  0.529422  kcal/mol 


Repeat  (y/n)?  y 


How  many  randomizations:  40 

What  fraction  of  the  distance  to  the  global  min  should  be  used  [0,1]:  0 
Use  how  many  points  from  previous  run  (global  only/all):  g 


Within  what  radius  of  phi[l] 
Within  what  radius  of  phi[2] 
Within  what  radius  of  phi[3] 
Within  what  radius  of  phi[4] 
Within  what  radius  of  phi[5] 
Within  what  radius  of  phi[6] 
Within  what  radius  of  phi[7] 
Within  what  radius  of  phi[8] 
Within  what  radius  of  phi{9] 


165.3140  (degs):  20 
177.6619  (degs):  5 
297.8685  (degs):  5 
295.7168  (degs):  30 
192.8867  (degs):  40 
295.6755  (degs):  30 
297.3322  (degs):  30 
177.1635  (degs):  30 
180.1444  (degs):  5 


— >  Using  40  randomizations 


running  QN-BFGS  (using  only  global  minimum) ...  pass  #4 


Elapsed  Real  Time  for  QN-BFGS  =  404.000  Seconds 
— >  42  points  reduced  to  42  unique  conformers 


solving  LP  ... 


Conformations:  42  (possibly  non-unique) 

Dihedral  Angles:  9 

Quad  Fit  Error  =  0.014822  kcal/mo!  (sum  of  errors) 
LP  Pivots  =  59 

Interpolated  Conformations: 

4  7  12  13  19  25  30  31  35  38  40  42  (12  out  of  42) 


Convex  Quadratic  Underestimator: 


c(0]  =  +325.080069 
c[l]  = +0.008313 
c[2]  =  -135.613593 
c[3]  =  +0.110036 
c[4]  =  +0.065432 
c[5]  = +0.01 1291 
c[6]  =  +0.040999 
c[7]  = +0.107349 
c[8]  =  -77.185222 
c[9]  =  -0.049789 


d[l]  =  +0.000000 
d[2]  = +43.781435 
d[3]  =  +0.000000 
d[4]  =  +0.000000 
d[5]  =  +0.000000 
d[6]  =  +0.000000 
d[7]  =  +0.000000 
d[8]  =  +24.898707 
d[9]  =  +0.000000 


Max  error  in  Quad  Fit  =  0.002218  kcal/mol 


Distribution  of  Errors  (%  of  Max) 


10%  20%  30%  40%  50%  60%  70% 

25  8  1  2  2  0  1 


Ranges  on  the  Dihedrals  (Degrees): 


179.779000  <=  phi[l]  <= 
177.013000  <=  phi[2]  <= 
297.070900  <=  phi[3]  <= 
295.240900  <=  phi[4]  <= 
169.495500  <=phi(5]  <= 
295.366900  <=  phi[6]  <= 
297.076300  <=  phi[7]  <= 
176.903600  <=  phi[8]  <= 
179.817100  <=  phi[9j  <= 


180.801700 

178.108200 

297.976400 

296.027200 

170.144400 

296.239900 

297.772100 

178.036600 

180.388500 


90%  100% 

0  2 


Elapsed  Real  Time  for  LP  =  1 .000  Seconds 


LP  Predicted  Dihedrals  Current  Global  Minimum 


phif  I  ]  =  3. 144677  rads*  (180. 176700  degs) 
phi[2]  =  3.097514  rads  (177.474460  degs) 
phi[3J  =  5.191054  rads*  (297.425500  degs) 
phi[4]  =  5.162241  rads*  (295.774600  degs) 
phi[5)  =  2.965339  rads*  (169.901400  degs) 


phif  I  ]  =  3. 144677  rads  (1 80. 1 76700  degs) 
phi[2]  =  3.096150  rads  (177.396300  degs) 
phi[3]  =  5.191054  rads  (297.425500  degs) 
phi[4]  =  5.162241  rads  (295.774600  degs) 
phi[5]  =  2.965339  rads  (169.901400  degs) 
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phi[6]  =  5.162111  rads*  (295.767200  degs)  phi[6]  =  5.162111  rads  (295.767200  degs) 

phi[7]  =  5.190196  rads*  (297.376300  degs)  phi[7]  =  5.190196  rads  (297.376300  degs) 

phi[8]  =  3.099969  rads  (177.615145  degs)  phi[8]  =  3.096898  rads  (177.439200  degs) 

phi(9]  =  3.144771  rads*  (180. 182 100  degs)  phi[9]  =  3  144771  rads  ( 1 80. 1 82 1  TO  degs) 

Energy  =  -3.007255  kcal/mol 

LP  Predic..on/Global  Minimum  Difference  (inf_norm)  =  0.003071  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  0.000158  kcal/mol 


Successive  LP  predicted  dihedrals  change  (inf-norm)  =  22.985300  degs 
Successive  global  conformers  change  (inf-norm)  =  0.275700  degs 
Energy  improvement  =  0  000983  kcal/mol 


Repeat  (y/n)?  y 

How  many  randomizations:  20 

What  fraction  of  the  distance  to  the  global  min  should  be  used  [0,1):  1 
Use  how  many  points  from  previous  run  (global  only/all):  g 

Within  what  radius  of  phi[l]  =  180.1767  (degs):  10 
Within  what  radius  of  phi[2]  =  177.3963  (degs):  10 
Within  what  radius  of  phi[3]  =  297.4255  (degs):  10 
Within  what  radius  of  phi[4]  =  295.7746  (degs):  10 
Within  what  radius  of  phi[5]  =  169.9014  (degs):  10 
Within  what  radius  of  phi[6]  =  295.7672  (degs):  10 
Within  what  radius  of  phi[7]  =  297.3763  (degs):  10 
Within  what  radius  of  phi[8]  =  177.4392  (degs):  10 
Within  what  radius  of  phi[9]  =  180.1821  (degs):  10 

-->  Using  20  randomizations 

running  QN-BFGS  (using  only  global  minimum) ...  pass  #5 
Elapsed  Real  Time  for  QN-BFGS  =  159.000  Seconds 
-->  22  points  reduced  to  22  unique  conformers 


solving  LP ... 


Conformations:  22  (possibly  non-unique) 

Dihedral  Angles:  9 

Quad  Fit  Error  =  0.008593  kcal/mol  (sum  of  errors) 
LP  Pivots  =  60 


Interpolated  Conformations: 
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1  34789  12  13  1920  (10outof22) 

Convex  Quadratic  Underestimator: 
c[0]  =  +0.000000 

c[l]  = +0.018242  d[l]  =  +0.000000 

c[2]  =  +0. 1 83987  d[2]  =  +0.000000 

c[3]  = +0.178514  d[3]  = +0.000000 

c[4]  =  +0.083497  d[4]  =  +0.000000 

c[5]  =  -0. 1 2472 1  d[5]  =  +0.000000 

c[6]  = -2.012263  d[6]  = +0.407210 

c[7]  =  +0. 1 3 1 898  d[7]  =  +0.000000 

c[8]  =  -0.060945  d[8]  =  +0.000000 

c[9]  =  -0.049604  d[9]  =  +0.000000 

Max  error  in  Quad  Fit  =  0.002296  kcal/mol 

Distribution  of  Errors  (%  of  Max) 

10%  20%  30%  40%  50%  60%  70%  80%  90%  100% 

13  4  0  1  2  0  0  1  1 

Ranges  on  the  Dihedrals  (Degrees): 

179.959300  <=  phi[l]  <=  180.312200 
177.251800  <=  phi [2]  <=  177.660500 
297.214300  <=  phi(3]  <=  297.691900 
295.571800  <=  phi[4]  <=  295.963200 
1 69.559600  <=  phi[5]  <=  170.075500 
295.394000  <=  phi[6]  <=  296.500900 
297.208700  <=  phi[7]  <=  297.696600 
177.086200  <=  phi{8]  <=  177.760500 
1 80.000700  <=  phi [9]  <=  1 80.53 1 600 


Elapsed  Real  Time  for  LP  =  0.000  Seconds 


LP  Predicted  Dihedrals  Current  Global  Minimum 

rhi[l)  =  3.144677  rads*  (180.176700 degs)  phi[l]  ^  3.144677  rads  (180.176700 degs) 

phi[2]  =  3.096150  rads*  (177.396300  degs)  phi[2]  =  3.096150  rads  (177.396300  degs) 

phi[3]  =  5.191054  rads*  (297.425500  degs)  phi[3J  =  5.191054  rads  (297.425500  degs) 

phi [4]  =  5.162241  rads*  (295.774600  degs)  phi[4)  =  5.162241  rads  (295.774600  degs) 

phi[5]  =  2.965339  rads*  (169.901400  degs)  phi(5]  =  2.965339  rads  (169.901400  degs) 

phi[6]  =  4.941585  rads  (283.131989  degs)  phi[6]  =  5.162111  rads  (295.767200  degs) 

phi[7]  =  5.190196  rads*  (297.376300  degs)  phi[7]  =  5.190196  rads  (297.376300  degs) 

phi[8]  =  3.096898  rads*  (177.439200  degs)  phi[8]  =  3.096898  rads  (177.439200  degs) 

phi[9)  =  3.144771  rads*  (180. 182 100  degs)  phi[9]  =  3.144771  rads  (180.182100  degs) 

Energy  =  -3.007255  kcal/mol 
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LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  0.220526  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  0.009902  kcal/moi 


Successive  LP  predicted  dihedrals  change  (inf-norm)  =  12.635200  degs 
Successive  global  conformers  change  (inf-norm)  =  0.000000  degs 
Energy  improvement  =  0.000000  kcal/mo! 


Repeat  (y/n)?  n 
Run  Complete.... 

E.  Global  Underestimator  Run:  22  beads 

Host  Machine:  csserverc  Wed  Apr  20  07: 12:32  EDT  1994 

How  many  randomizations  (at  least  39):  100 

Within  what  radius  of  phifl]  =  330.0000  (degs):  180 
Within  what  radius  of  phi[2]  =  278.0000  (degs):  180 
Within  what  radius  of  phi[3]  =  123.0000  (degs):  180 
Within  what  radius  of  phi[4]  =  31.0000  (degs):  180 
Within  what  radius  of  phi[5]  =  160.0000  (degs):  180 
Within  what  radius  of  phi[6]  =  97.0000  (degs):  180 
Within  what  radius  of  phi[7]  =  27.0000  (degs):  180 
Within  what  radius  of  phi[8]  =  232.0000  (degs):  180 
Within  what  radius  of  phi[9]  =  38.0000  (degs):  180 
Within  what  radius  of  phi[10]  =  213.0000  (degs):  180 
Within  what  radius  of  plu[ll]  =  311.0000  (degs):  180 
Within  what  radius  of  phi[12]  =  137.0000  (degs):  180 
Within  what  radius  of  phi[  13]  =  28.0000  (degs):  180 
Within  what  radius  of  phi[  14]  =  16.0000  (degs):  180 
Within  what  radius  of  phi[  15]  =  323.0000  (degs):  180 
Within  what  radius  of  phi[16]  =  222.0000  (degs):  180 
Within  what  radius  of  phi[17]  =  289.0000  (degs):  180 
Within  what  radius  of  phi[  18]  =  120.0000  (degs):  180 
Within  what  radius  of  phi[  19]  =  179.0000  (degs):  180 

-->  Using  100  randomizations 

running  QN-BFGS  (using  all  local  minima) ...  pass  #1 

Elapsed  Real  Time  for  QN-BFGS  =  14465.000  Seconds 

— >  99  points  reduced  to  99  unique  conformers 


solving  LP  ... 
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Conformations:  99  (possibly  non-unique) 

Dihedral  Angles:  19 

Quad  Fit  Error  =  261.999901  kcal/mol  (sum  of  errors) 

LP  Pivots  =  126 

Interpolated  Conformations: 

2  4  5  7  17  1825  26  28  30  33  35  38  39  4042  45  52  53  56  57  59  6063  64  71  72  76  77  7981  83  84 
87  94  99  (36  out  of  99) 


Convex  Quadratic  Underestimator: 


c[0]  = +15.809699 
c[l]  =  -0.238399 
c[2]  = -1.555963 
c[3]  = +0.041718 
c[4]  =  +0.320074 
c[5]  =  -0.523003 
c[6]  = -1.399913 
c[7]  =  -2.284349 
c[8]  =  -0.884438 
c[9]  =  -0.529013 
c[10]  = -1.135499 
c[ll]  =  +0.393577 
c[12]  = -1.486392 
c[13]  = +0.095308 
c[14]  = -0.221887 
c[15]  =  -2.938166 
c[16]  = -2.479375 
c[17]  = -1.420245 
c[  18]  =-0.300846 
c[19]  = -1.454993 


d[l]  = +0.092840 
d[2]  = +0.3741 13 
d[3]  =  +0.000000 
d[4]  = +0.148104 
d[5]  = +0.140036 
d[6]  =  +0.504272 
d[7]  =  +0.802268 
d[8]  =  +0.334472 
d[9]  =  +0.000000 
d[10]  =  +0.241783 
d[ll]  = +0.000000 
d[12]  =  +0.398304 
d[13]  = +0.014153 
d[  14]  =  +0.082623 
df  15]  =  +0.789750 
d[16]  = +0.895932 
d[17]  =  +0.501218 
d[18]  =  +0.192963 
d[19]  =  +0.375462 


Max  error  in  Quad  Fit  =  10.712270  kcal/mol 


Distribution  of  Errors  (%  of  Max) 


10%  20%  30%  40%  50%  60%  70% 

43  10  8  12  9  4  6 


80%  90%  100% 

3  2  2 


Ranges  on  the  Dihedrals  (Degrees): 


60.764300 

58.390200 

55.255200 

57.929100 

54.104100 

61.046500 

53.912900 

51.561400 

53.754300 

51.790700 

52.621200 

52.233600 


<=  phi[l]  <=  300.277800 
<=  phi[2]  <=  301.986700 
<=  phi[3]  <=  300.553900 
<=  phi[4]  <=  302.100000 
<=phi[5]<=  302.815000 
<=  phi[6]  <=  304.852100 
<=phi[7]<=  309.291800 
<=  phi[8]  <=  305.891100 
<=  phi(9]  <=  307.740200 
<=  phi]  10]  <=  307.421300 
<=  phi]  1 1  ]  <=  310.690700 
<=  phi]  12]  <=  305.259600 
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59.822800  <=  phi[13]  <=  306.844000 
54.681400  <=  phi[  14]  <=  303.379200 
51.037900  <=  phi[15]  <=  303.680100 
57.200900  <=  phi[16]  <=  304.314400 
57.143800  <=  phi [17]  <=  301.433200 
60.3 1 5 1 00  <=  phi[  1 8]  <=  300.888600 
59.358800  <=  phi[19]  <=  300.052300 


Elapsed  Real  Time  for  LP  =  4.000  Seconds 


LP  Predicted  Dihedrals  Current  Global  Minimum 


phi[l]  =  2.567848  rads  (147.126848  degs) 
phi[2]  =  4.159072  rads  (238.297287  degs) 
phi[3]  =  1.078123  rads*  (61.771900  degs) 
phi[4]  =  4. 122042  rads  (236. 1 75597  degs) 
phi[5]  =  3.734775  rads  (213.986865  degs) 
phi[6]  =  2.776107  rads  (159.05921 1  degs) 
phi[7]  =  2.847364  rads  (163.141938  degs) 
phi[8]  =  2.644281  rads  (151.506149  degs) 
phi[9]  =  5.148346  rads*  (294.978500  degs) 
phi[  10]  =  4.696356  rads  (269.081368  degs) 
phi[  11]  =  5.165801  rads*  (295.978600  degs) 
phi[  1 2]  =  3.731803  rads  (213.816553  degs) 
phi[13]  =  5.832251  rads  (334.163347  degs) 
phi[  14]  =  2.685536  rads  (153.869850  degs) 
phi[  15]  =  3.720375  rads  (213.161774  degs) 
phi[16]  =  2.767370  rads  (158.558600  degs) 
phi[17]  =  2.833587  rads  (162.352598  degs) 
phil  1 8]  =  1 .559086  rads  (  89.329074  degs) 
phi[19]  =  3.875207  rads  (222.033010  degs) 


phi[l]  =  3.105243  rads  (177.917300  degs) 
phi[2]  =  3.183406  rads  (182.395700  degs) 
phi[3]  =  1.078123  rads  (  61.771900  degs) 
phi[4]  =  3. 164792  rads  (181 .329200  degs) 
phi[5)  =  3.149208  rads  (180.436300  degs) 
phi[6]  =  1.088523  rads  (  62.367800  degs) 
phi[7]  =  3.071958  rads  (176.010200  degs) 
phi[8]  =  3.145982  rads  (180.251500  degs) 
phi[9]  =  5.148346  rads  (294.978500  degs) 
phi[  10]  =  3.022744  rads  (173.190500  degs) 
phi[  11]  =  5.165801  rads  (295.978600  degs) 
phi[  12]  =  5. 196676  rads  (297.747600  degs) 
phi[13]  =  5.129713  rads  (293.910900  degs) 
phi[  14]  =  1.678839  rads  ( 96.190400  degs) 
phi[  15]  =  1 .144625  rads  ( 65.582200  degs) 
phi[16]  =  3.288972  rads  (188.444200  degs) 
phi[17]  =  1.174315  rads  (  67.283300  degs) 
phi[18]  =  3.072954  rads  (176.067300  degs) 
phi[  19]  =  3.208947  rads  (183.859100  degs) 

Energy  =  -7.449145  kcal/mol 


LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  2.575749  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  4.5954 1 6  kcal/mol 


Repeat  (y/n)?  y 


How  many  randomizations:  75 

What  fraction  of  the  distance  to  the  global  min  should  be  used  [0,1]:  0 
Use  how  many  points  from  previous  run  (global  only/all):  g 


Within  what  radius 
Within  what  radius 
Within  what  radius 
Within  what  radius 
Within  what  radius 
Within  what  radius 
Within  what  radius 


of  phi[l]  =  147.1268  (degs):  35 
of  phi[2]  =  238.2973  (degs):  55 
of  phi[3]  =  61.7719  (degs):  180 
of  phi[4]  =  236. 1756  (degs):  55 
of  phi[5]  =  213.9869  (degs):  35 
of  phi[6]=  159.0592  (degs):  100 
of  phi[71  =  163.1419  (degs):  15 


85 


Within  what  radius  of  phi[8]  =  151.5061  (degs):  30 
Within  what  radius  of  phi[9]  =  294.9785  (degs):  180 
Within  what  radius  of  phi[  10]  =  269.0814  (degs):  100 
Within  what  radius  of  phi[  11]  =  295.9786  (degs):  180 
Within  what  radius  of  phi [12]  =  213.8166  (degs):  85 
Within  what  radius  of  phi[13]  =  334.1633  (degs):  35 
Within  what  radius  of  phi[  14]  =  153.8699  (degs):  60 
Within  what  radius  of  phi[  15]  =  213.1618  (degs):  150 
Within  what  radius  of  phi[16]  =  158.5586  (degs):  30 
Within  what  radius  of  phi[17]  =  162.3526  (degs):  105 
Within  what  radius  of  phi[18]  =  89.3291  (degs):  90 
Within  what  radius  of  phi[  19]  =  222.0330  (degs):  40 

— >  Using  75  randomizations 

running  QN-BFGS  (using  only  global  minimum) ...  pass  #2 


Elapsed  Real  Time  for  QN-BFGS  =  6746.000  Seconds 
-->  77  points  reduced  to  77  unique  conformers 


solving  LP ... 


Conformations:  77  (possibly  non-unique) 

Dihedral  Angles:  19 

Quad  Fit  Error =180.596223  kcal/mol  (sum  of  errors) 
LP  Pivots  =  1 59 


Interpolated  Conformations: 

1  6  8  9  21  23  24  26  30  31  34  41  42  44  45  48  51  60  62  64  67  69  72  75  76  (25  out  of  77) 


Convex  Quadratic  Underestimator: 


c[0]  =  +144.561300 
c[l]  = -34.584215 
c[2]  = -13.561606 
c[3]  =  -0.939967 
c[4]  =  +0.845433 
c[5]  =  -31.655626 
c[6]  =  -0.457572 
c[7]  = -5.701583 
c[8]  = -1.826113 
c(9]  =  -0.438230 
c[  10]  =  +0.490778 
c[  11]  =  -2.679081 
c[12]  =  -0.166776 
c[13]  = -0.366454 
c(14]  = +0.147614 
c[15]  = +0.104883 


d[l]  = +16.810994 
d[2]  =  +3.534313 
d[3]  =  +0.000000 
d[4]  =  +0.000000 
d[5]  =  +7.783377 
d[6]  =  +0.000000 
d[7]  =  +0.000000 
d[8]  =  +0.000000 
d{9]  = +0.000000 
d[  10]  =  +0.000000 
d[ll]  =  +0.807040 
d[12]  = +0.000000 
d[  1 3]  = +0.000000 
d[  14]  =  +0.000000 
d[  1 5]  =  +0.000000 
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c[16]  =  -3.023704 
c[17]  = -0.722325 
c[18]  =  -0.905696 
cfl9]  = +0.756502 


d[  16]  =  +0.000000 
d[17]  = +0.227437 
d[18]  = +0.000000 
d[ !  9}  =  +0.000000 


Max  error  in  Quad  Fit  =  9.905244  kcal/mol 

Distribution  of  Errors  (%  of  Max) 

10%  20%  30%  40%  50%  60%  70%  80%  90%  100% 

30  6  10  12  7  7  3  1  0  1 

Ranges  on  the  Dihedrals  (Degrees): 

61.225300  <=  phi[  1  ]  <=  181.832900 
175.450700  <=phi[2]  <=  301.786100 
60.548900  <=  phi[3]  <=  301.842700 
169.540400  <=  phi [4]  <=  299.865600 
159.878900  <=  phi]5]  <=  299.666000 
53.564400  <=  phi(6]  <=  303.470100 
165.506600  <=  phi[7]  <=  207.251 100 
75.822700  <=  phi[8]  <=  205.201600 
57.909800  <=  phi[9]  <=  307. 113100 
55.03 1 500  <=  phi[  1 0]  <=  3 1 1 .93 1 900 
49.647200  <=  phi]  1 1  ]  <=  301 .752600 
1 50.4 1 3300  <=  phi[  1 2]  <=  3 1 1 .839200 
55.901 100  <=  phi[13]  <=  308.383200 
58.148900  <=  phi[14]  <=  198.021200 
55.392300  <=  phi[15]  <=  305.325100 
170.785700  <=  phi[  16]  <=  193.669800 
57.015500  <=  phi[17]  <=  302.169300 
59.621500  <=  phi]  18]  <=  183.323300 
174.983100  <=  phi]  19]  <=  297.404000 


Elapsed  Real  Time  for  LP  =  3.000  Seconds 


LP  Predicted  Dihedrals  Current  Global  Minimum 


phi]  1  ]  =  2.057238  rads  (117.871053  degs) 
phi[2]  =  3.837126  rads  (219.851 153  degs) 
phi]3]  =  1.056778  rads*  (  60.548900  degs) 
phi]4]  =  3.172085  rads*  (181.747100  degs) 
phi[5]  =  4.067081  rads  (233.026586  degs) 
phi{6]  =  5.234387  rads*  (299.908300  degs) 
phi]7]  =  3.154526  rads*  (180.741000  degs) 
phi]8]  =  3.083152  rads*  (176.651600  degs) 
phi]9]  =  5.186810  rads*  (297.182300  degs) 
phi]  10]  =  5.233930  rads*  (299.882100  degs) 
phi]  1 1  ]  =  3.319638  rads  (190.201272  degs) 
phi]  12]  =  5.244374  rads*  (300.480500  degs) 
phi]  13]  =  5.21 1585  rads*  (298.601800  degs) 


phi[l]  =  1.101978  rads  (  63.138700  degs) 
phi[2]  =  3. 1 576 1 8  rads  ( 1 80.9 1 8200  degs) 
phi]  3]  =  1 .056778  rads  (  60.548900  degs) 
phi[4]  =  3.172085  rads  (181.747100  degs) 
phi[5]  =  3.115551  rads  (178.507900  degs) 
phi[6]  =  5.234387  rads  (299.908300  degs) 
phi]7]  =  3.154526  rads  (180.741000  degs) 
phi [8]  =  3.083152  rads  (176.651600  degs) 
phi[9]  =  5.186810  rads  (297. 182300  degs) 
phi]  10]  =  5.233930  rads  (299.882100  degs) 
phi]  1 1  ]  =  3.052165  rads  (174.876200  degs) 
phi]  12]  =  5.244374  rads  (300.480500  degs) 
phi]  13]  =  5.21 1585  rads  (298.601800  degs) 
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phi[  14]  =  3.067434  rads*  ( 175.751 000  degs) 
phi[15]  =  5.241103  rads*  (300.293100  degs) 
phi[16]  =  3.108295  rads*  (178.092200  degs) 
phif  1 7]  =  3. 1 75934  rads  (181 .967639  degs) 
phi[18]  =  1.040591  rads*  (  59.621500  degs) 
phif  19]  =  3.149848  rads*  (180.473000  degs) 


phi[14]  =  3.067434  rads  (175.751000  degs) 
phif  15]  =  5.241103  rads  (300.293100  degs) 
phif  16]  =  3.108295  rads  (178.092200  degs) 
phif  1 7]  =  3. 1 55606  rads  ( 1 80.802900  degs) 
phif  18]  =  1.040591  rads  (  59.621500  degs) 
phif  19]  =  3.149848  rads  (180.473000  degs) 


Energy  =  -8.726705  kcal/mol 


LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  0.955260  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  12.038637  kcal/mol 


Successive  LP  predicted  dihedrals  change  (inf-norm)  =  140.849100  degs 
Successive  global  conformers  change  (inf-norm)  =  237.540500  degs 
Energy  improvement  =  1 .277560  kcal/mol 


Repeat  (y/n)?  y 

How  many  randomizations:  65 

What  fraction  of  the  distance  to  the  global  min  should  be  used  [0,1]:  0.5 
Use  how  many  points  from  previous  run  (global  only/all):  g 

Within  what  radius  of  phif  1  ]  =  90.5049  (degs):  30 
Within  what  radius  of  phi[2]  =  200.3847  (degs):  20 
Within  what  radius  of  phi[3]  =  60.5489  (degs):  180 
Within  what  radius  of  phi[4]  =  181.7471  (degs):  60 
Within  what  radius  of  phif5]  =  205.7672  (degs):  25 
Within  what  radius  of  phi[6]  =  299.9083  (degs):  60 
Within  what  radius  of  phi[7]  =  180.7410  (degs):  5 
Within  what  radius  of  phi[8]  =  176.6516  (degs):  10 
Within  what  radius  of  phi[9]  =  297.1823  (degs):  60 
Within  what  radius  of  phif  10]  =  299.8821  (degs):  10 
Within  what  radius  of  phif  1 1  ]  =  182.5387  (degs):  5 
Within  what  radius  of  phif  12]  =  300.4805  (degs):  60 
Within  what  radius  of  phif  13]  =  298.6018  (degs):  5 
Within  what  radius  of  phif  14]  =  175.7530  (degs):  60 
Within  what  radius  of  phif  15]  =  300.2931  (degs):  60 
Within  what  radius  of  phif  16]  =  178.0922  (degs):  10 
Within  what  radius  of  phi[17]  =  181.3853  (degs):  5 
Within  what  radius  of  phif  1 8]  =  59.62 1 5  (degs):  60 
Within  what  radius  of  phif  19]  =  180.4730  (degs):  20 

— >  Using  65  randomizations 

running  QN-BFGS  (using  only  global  minimum) ...  pass  #3 


Elapsed  Real  Time  for  QN-BFGS  =  7484.000  Seconds 


— >  67  points  reduced  to  67  unique  conformers 


solving  LP ... 


Conformations:  67  (possibly  non-unique) 

Dihedral  Angles:  19 

Quad  Fit  Error  =  0.713424  kcal/mol  (sum  of  errors) 
LP  Pivots  =  124 


Interpolated  Conformations: 


1  2  3  6  8  12  17  21  24  26  29  30  39  40  41  42  43  44  46  47  48  49  50  52  55  58  61  66  67  (29  out  of  67) 
Convex  Quadratic  Underestimator: 


c[0]  =  +862.353679 
c[l]  = +1.576256 
c[2]  =  +11.505862 
c[3]  =  -0.755355 
c[4]  =  +0.861950 
c[5]  = +1.905837 
c[6]  =  -14.157930 
c[7]  =  -8.344343 
c[8]  =  -3.335837 
c[9]  =  +3.106748 
c[10]  = -68.281811 
c[ll]  =  -8.726862 
c[12]  = -3.075571 
c[13]  = +3.623808 
c[14]  = -22.205927 
c[15]  = +0.243679 
c[16]  = -31.923569 
c[17]  = -0.710224 
c[18]  =  -5.924973 
c[19]  =  -391.635793 


d[l]  = +0.000000 
d[2]  =  +0.000000 
d[3]  =  +0.265304 
d[4]  =  +0.000000 
d[5]  =  +0.000000 
d[6]  =  +4.433968 
d[7]  =  +0.386241 
d[8]  =  +0.000000 
d[9]  =  +0.000000 
d[10]  =  +13.059085 
dfll]  =  +2.753363 
d[  1 2]  = +0.000000 
d[  1 3]  = +0.000000 
d[14]  =  +7.411516 
d(  15]  = +0.000000 
d[16]  = +10.592633 
d[17]  =  +0.000000 
d[18]  =  +1.949442 
d[19]  = +126.049962 


Max  error  in  Quad  Fit  =  0.073463  kcal/mol 

Distribution  of  Errors  (%  of  Max) 

10%  20%  30%  40%  50%  60%  70%  80%  90%  100% 

41  6  6  6  5  0  1  0  1  1 

Ranges  on  the  Dihedrals  (Degrees): 

59.668300  <=  phi[l]  <=  180.233200 
173.662000  <=  phi[2]  <=  186.797300 
60.465200  <=  phi[3]  <=  302.153700 
168.498000  <=  phi[4]  <=  302.173200 
153.880200  <=  phi[5]  <=  297.450400 
58.388000  <=  phi(6]  <=  306.578200 


173.159200  <= 
64.552300  <= 
72.285500  <= 
1 88.202400  <= 
161.972500  <= 
177.360700  <= 
57.296400  <= 
60.468000  <= 
176.636700  <= 
64.959300  <= 
173.914200  <= 
55.501000  <= 
176.297300  <= 


phi[7]  <= 
phi[8]  <= 
phi  [9]  <= 
phi[  10]  < 
phi[ll]  < 
phi[12]  < 
phi[  1 3]  <; 
phi[  14]  <; 
phif  15]  <: 
phi[16]  < 
phif  17]  c 
phif  18]  c 
phif  19] « 


302.506900 
294.232600 
303.774700 
=  305.621100 
=  294.907100 
=  306.147500 
301.846000 
293.192800 
:=  304.215300 
301.386800 
:=  300.148500 
302.929000 
:=  184.395100 


Elapsed  Real  Time  for  LP  =  3.000  Seconds 


LP  Predicted  Dihedrals  Current  Global  Minimum 


phif  1  ]  =  1.114948  rads*  (  63.881800  degs) 
phi[2]  =  3.146228  rads*  (180.265600  degs) 
phif3]  =  2.847130  rads  (163.128538  degs) 
phi[4]  =  3.074986  rads*  (176.183700  degs) 
phi[5]  =  5.191489  rads*  (297.450400  degs) 
phi[6]  =  3.193061  rads  (182.948915  degs) 
phi[7]  =  2.754424  rads  (157.816899  degs) 
phif8]  =  5.135328  rads*  (294.232600  degs) 
phi [9]  =  5. 124215  rads*  (293.595900  degs) 
phiflO]  =  5.228683  rads  (299.581448 degs) 
phif  1 1  ]  =  3 . 1 69528  rads  ( 1 8 1 .600596  degs) 
phif  12]  =  5.130233  rads*  (293.940700  degs) 
phif  13]  =  5.268207  rads*  (301.846000  degs) 
phif  14]  =  2.996138  rads  (171.666080  degs) 
phif  15]  =  5.266020  rads*  (301.720700  degs) 
phi[16]  =  3.013752  rads  (172.675271  degs) 
phif  17]  =  3.223166  rads*  (184.673800  degs) 
phif  18]  =  3.039320  rads  (174.140207  degs) 
phif  19]  =  3.106989  rads  (178.017333  degs) 


phif  1]  =  1.114948  rads  (  63.881800  degs) 
phi[2]  =  3. 146228  rads  ( 1 80.265600  degs) 
phi[3]  =  5.212571  rads  (298.658300  degs) 
phi[4]  =  3.074986  rads  (176.183700  degs) 
phi[5]  =  5.191489  rads  (297.450400  degs) 
phi[6]  =  5.151067  rads  (295.134400  degs) 
phi[7]  =  3.022198  rads  (173.159200  degs) 
phif 8]  =  5.135328  rads  (294.232600  degs) 
phi[9]  =  5.124215  rads  (293.595900  degs) 
phiflO]  =  5.144740  rads  (294.771900  degs) 
phif  1 1]  =  2.983309  rads  (170.931000  degs) 
phif  12]  =  5.130233  rads  (293.940700  degs) 
phif  13]  =  5.268207  rads  (301.846000  degs) 
phif  1 4]  =  3. 1 72427  rads  (181 .766700  degs) 
phif  1 5]  =  5.266020  rads  (301 .720700  degs) 
phif  16]  =  3.190450  rads  (182.799300  degs) 
phif  17]  =  3.223166  rads  (184.673800  degs) 
phif  18]  =  1.153593  rads  (  66.096000  degs) 
phif  1 9]  =  3. 1 69434  rads  (181 .595200  degs) 


Energy  =  -10.974719  kcal/mol 


LP  Prediction/Global  Minimum  Difference  (inf_norm)  =  2.365441  degs 
LP  Prediction/Global  Minimum  Quad  Function  Difference  =  11.392121  kcal/mol 


Successive  LP  predicted  dihedrals  change  (inf-norm)  =  1 17.581000  degs 
Successive  global  c-mformers  change  (int-norm)  =  238.109400  degs 
Energy  improvement  =  2.248014  kcal/mol 


Repeat  (y/n)?  y 


How  many  randomizations:  60 

What  fraction  of  the  distance  to  the  global  min  should  be  used  [0,1  ]:  0 
Use  how  many  points  from  previous  run  (global  only/all):  g 


Within  what  radius  of  phi[l]  = 
Within  what  radius  of  phi  [2]  = 
Within  what  radius  of  phi[3]  = 
Within  what  radius  of  phi[4]  = 
Within  what  radius  of  phi[5]  = 
Within  what  radius  of  phi[6]  = 
Within  what  radius  of  phi[7]  = 
Within  what  radius  of  phi[8]  = 
Within  what  radius  of  phi[9]  = 
Within  what  radius  of  phi[  10] 
Within  what  radius  of  phi[  11] 
Within  what  radius  of  phi[12] 
Within  what  radius  of  phi[  13] 
Within  what  radius  of  phi[14] 
Within  what  radius  of  phi[  15] 
Within  what  radius  of  phi[  16] 
Within  what  radius  of  phi[  17] 
Within  what  radius  of  phi[18] 
Within  what  radius  of  phi[  19] 

~>  Using  60  randomizations 


63.8818  (degs):  5 
180.2656  (degs):  5 
163.1285  (degs):  100 
176.1837  (degs):  5 
297.4504  (degs):  45 
182.9489  (degs):  100 
157.8169  (degs):  15 
294.2326  (degs):  5 
293.5959  (degs):  20 
=  299.5814  (degs):  5 
=  181.6006  (degs):  10 
=  293.9407  (degs):  20 
=  301.8460  (degs):  5 
=  171.6661  (degs):  10 
=  301.7207  (degs):  10 
=  172.6753  (degs):  10 
=  184.6738  (degs):  5 
=  174.1402  (degs):  100 
=  178.0173  (degs):  5 


running  QN-BFGS  (using  only  global  minimum) ...  pass  #4 


Elapsed  Real  Time  for  QN-BFGS  =  5461.000  Seconds 
~>  62  points  reduced  to  62  unique  conformers 


solving  LP ... 


Conformations:  62  (possibly  non-unique) 

Dihedral  Angles:  19 

Quad  Fit  Error  =4.658270  kcal/mol  (sum  of  errors) 

LP  Pivots  =  179 

Interpolated  Conformations. 

1  69  12  15  17  21  22  24  30  33  36  38  41  48  49  50  51  55  56  57  59  60  (23  out  of  62) 
Convex  Quadratic  Underestimator: 


c[0]  = +1555.161908 
c[l]  =  -34.543558 


d[l]  =  +0.000000 


c[2]  =  -223.829278 
c[3]  = -0.293176 
C[4]  =  -4.439108 
c[5]  =  +4.849708 
c[6]  = -2.616177 
c[7]  = +13.611229 
c[8]  =  -6.794006 
c[9]  = -14.933655 
c[  10]  =  +2.805539 
c[ll]  = -637.240473 
c[12]  = -19.966918 
c[13]  = -1.045903 
c[14]  =  -121.797891 
c[15]  =  -13.684770 
c[16]  =  +8.283929 
c[17]  = +26.833230 
c[18]  = +0.768751 
c[19]  = +30.518616 


d[2]  =  +75.446077 
d[3]  =  +0.000000 
d[4]  =  +0.000000 
d[5]  =  +0.000000 
d[6]  =  +0.000000 
d[7]  =  +0.000000 
d[8]  =  +0.000000 
d[9]  =  +0.000000 
d[  10]  = +0.000000 
d[ll]  =  +210.926966 
d[  1 2]  = +0.000000 
d[  1 3]  = +0.000000 
d[  1 4]  =  +35.933857 
d[15]  =  +0.000000 
d[  16]  = +0.000000 
d[  1 7]  = +0.000000 
d[  1 8]  = +0.000000 
d[  1 9]  = +0.000000 


Max  error  in  Quad  Fit  =  0.462053  kcal/mol 


Distribution  of  Errors  (%  of  Max) 
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Ranges  on  the  Dihedrals  (Degrees): 


59.754100  <= 
174.363600  <= 
58.142400  <= 
165.275500  <= 
1 83.900100  <= 
52.505800  <= 
168.653900  <= 
288.604600  <= 
290.888000  <= 
289.843600  <= 
170.113300  <= 
289.565600  <= 
291.950200  <= 
168.734700  <= 
291.373000  <= 
172.934100  <= 
173.068600  <= 
57.784600  <= 
175.891700  <= 


phi[  1  ]  <=  64.670100 
:  phi[2]  <=  185.634700 
phi[3]  <=  299.767500 
:  phi[4]  <=  193.215700 
:phi[5]<=  301.015500 
phi[6]  <=  303.772100 
phi[7]  <=  187.344800 
phi[8]  <=  303.332400 
phi[9]  <=  300.782500 
phi]  10]  <=  306.124200 
phi]  11]  <=  182.814400 
phi]  12]  <=  306.767000 
phi]  13]  <=  303.437500 
phi]  14]  <=  192.165600 
phi]  15]  <=  307.080400 
phi]  16]  <=  186.826300 
phi]  17]  <=  191.041000 
phi]  18]  <=  300.637400 
:  phi]  19]  <=  182.798600 
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Elapsed  Real  Time  for  LP  =  4.000  Seconds 


100% 

1 


LP  Predicted  Dihedrals 


Current  Global  Minimum 
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phi[  1]  =  1.1 14948  rads*  ( 63.881 800  degs) 
phi[2]  =  2.966745  rads  (169.981972  degs) 
phi[3]  =  5.212571  rads*  (298.658300  degs) 
phi{4]  =  3.074986  rads*  (176.183700  degs) 
phi[5]  =  5.191489  rads*  (297.450400  degs) 
phi[6]  =  5.151067  rads*  (295.134400  degs) 
phi[7]  =  3.022198  rads*  (173.159200  degs) 
phi[8]  =  5.135328  rads*  (294.232600  degs) 
phif9]  =  5.124215  rads*  (293.595900  degs) 
phi(10j  =  5.144740  rads*  (294.771900  degs) 
phi[  1 1  ]  =  3.02 1 1 43  rads  ( 1 73.098729  degs) 
phi(12]  =  5.130233  rads*  (293.940700  degs) 
phi[  13]  =  5.268207  rads*  (30!  .846000  degs) 
phi]  14]  =  3.389502  rads  (194.204177  degs) 
phi[15]  =  5.266020  rads*  (301.720700  degs) 
phi]  16]  =  3.190450  rads*  (182.799300  degs) 
phi]  1 7]  =  3.223 1 66  rads*  ( 1 84.673800  degs) 
phi]  1 8]  =  1 . 1 53593  rads*  (  66.096000  degs) 
phi]  19]  =  3.169434  rads*  (181.595200  degs) 


phi]  1  ]  =  1.114948  rads  (  63.881800  degs) 
phi[2]  =  3.146228  rads  (180.265600  degs) 
phi[3]  =  5.212571  rads  (298.658300  degs) 
phi[4]  =  3.074986  rads  ( 1 76. 1 83700  degs) 
phi[5]  =  5.191489  rads  (297.450400  degs) 
phi[6]  =  5.151067  rads  (295.134400  degs) 
phi[7]  =  3.022198  rads  (173.159200  degs) 
phi[8]  =  5.135328  rads  (294.232600  degs) 
phi[9]  =  5.124215  rads  (293.595900  degs) 
phi]  10]  =  5.144740  rads  (294.77 1900  degs) 
phi]  1 1  ]  =  2.983309  rads  (170.93 1000  degs) 
phi]  12]  =  5.130233  rads  (293.940700  degs) 
phi]  13]  =  5.268207  rads  (301. 846000  degs) 
phi]  14]  =  3.172427  rads  (181.766700  degs) 
phi]  15]  =  5.266020  rads  (301.720700  degs) 
phi]  1 6]  =  3. 1 90450  rads  ( 1 82.799300  degs) 
phi]  17]  =  3.223166  rads  (184.673800  degs) 
phi]  18]  =  1.153593  rads  (  66.096000  degs) 
phi]  1 9]  =  3. 1 69434  rads  (181 .595200  degs) 


Energy  =  -10.974719  kcal/mol 

LP  Prediction/Globa!  Minimum  Difference  (inf_norm)  =  0.217075  degs 
LP  Prcdiction/Global  Minimum  Quad  Function  Difference  =  2.212808  kcal/mol 


Successive  LP  predicted  dihedrals  change  (inf-norm)  =  135.529800  degs 
Successive  global  conformers  change  (inf-norm)  =  0.000000  degs 
Energy  improvement  =  0.000000  kcal/mol 

Repeat  (y/n)?  n 
Run  Complete... 

F.  Q'  as  a  Negative  Definite  Matrix 

As  stated  in  Chapter  VIII  (Section  C)  if  Q'  =  Q-\il,  where  p  =  I  +  ^max  and 

X.max  is  the  maximum  eigenvalue  of  Q,  then  Q '  is  negative  definite.  This  can  be  shown  as 
follows:  Let  k  and  v>  be  any  eigenvalue  /  eigenvector  pair  for  Q,  i.e. 

Qv  =  Xv  Equation  12.1 


Since  Q'  =  (?-p/,  then 


Q'v  = 
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=  (fi-d+^DP  =  (Qv-  (l  +  ^max)lv) 


Substituting  Equation  12.1, 

Q'v  =  [Xv- ( 1  +  ^max)  Ivj  =  Zv-v-\maxv 

=  |^-^majtJv-vJ  =  [*-*««- *]* 

since  X  <  \max ,  (  X-  \max  —  1 J  <  —1 .  Thus  all  eigenvectors  of  Q  are  also  eigenvectors  of 
Q 'and  all  eigenvalues  of  Q' are  negative.  Hence,  the  matrix  Q' is  negative  definite.  □ 


£Mt> 
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